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ABSTRACT 

An exponential finite difference algorithm, as first presented by 
Bhattacharya for one-dimensional unsteady state, heat conduction In 
Cartesian coordinates, has been extended. The finite difference 
algorithm developed was used to solve the diffusion equation In 
one-dimensional cylindrical coordinates and applied to two- and 
three-dimensional problems In Cartesian coordinates. The method was 
also used to solve nonlinear partial differential equations In one 
(Burger's equation) and two (Boundary Layer equations) dimensional 
Cartesian coordinates. Predicted results were compared to exact 
solutions where available, or to results obtained by other numerical 
methods. It was found that the exponential finite difference method 
produced results that were more accurate than those obtained by other 
numerical methods, especially during the Initial transient portion of 
the solution. Other applications made using the exponential finite 
difference technique Included unsteady one-dimensional heat transfer 
with temperature varying thermal conductivity and the development of 
the temperature field In a laminar Couette flow. 
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NOMENCLATURE 


,b^ ,c 1 

Thomas algorithm variables 

B 

Blot number 

C P 

material specific heat, J/kg • K (Btu/lb m • °F) 

h 

convection heat transfer coefficient, W/M 2 • °C 
(Btu/ft 2 • hr • °F ) 

I.J.k 

nodal location In x,y, and z spatial coordinate directions 
respectively 

Jo» J i 

Bessel functions of zero and first order 

k 

thermal conductivity, W/M • °C (Btu/hr • °F • ft) 

k? 

thermal conductivity at 1^ position, n^ time step, 

W/M • °C (Btu/hr • °F • ft) 

% 

L 

distance between plates, M (ft) 

M 

dimensionless drive number 

m 

number of sub-intervals 

N 

number of nodes In a spatial direction 

n 

time step position designation 

Q 

heat flux, W/M 2 (Btu/hr • ft 2 ) 

r 

spatial coordinate; cylindrical coordinates, M (ft) 

T 

temperature, °C (°F) 

t 

time, sec 

At 

time between time steps n and n + 1 

U 

Couette flow velocity. M/s (ft/s) 

x.y.z 

spatial coordinates, Cartesian coordinates, M (ft) 

AX.Ay.AZ 

distance between nodal positions In the x,y, and z spatial 
directions 

a 

thermal dlffuslvlty, M 2 /s (ft 2 /s) 

0 

rate of thermal conductivity variation 

Y 

At/pC p ( Ax) 2 ; (W/M • “C)' 1 ((Btu/hr • °F • ft 2 )' 1 ) 


v 

PRECEDING PM* 51 DLA? ' K NOT FILMED 


it constant 

constant used In exponential finite difference method with 
temperature varying thermal conductivity 

< x finite difference operator 

Xn, mth eigenvalue of Bessel function 

X?,v? Thomas algorithm variables dependent on time step and spatial 
location 

l amplification factor 

v kinematic viscosity, M 2 /s (ft 2 /s) 

p material density, kg/M 3 (lbm/ft 3 ) 

<p,i)»,e separation variables 

K - dimensionless time 

(Ax)' 


vl 


I. INTRODUCTION 


Partial differential equations have many Important applications In 
the fields of engineering and physics. Many exact solutions exist 
depending on, the partial differential equation, the boundary 
conditions, and the number of spatial dimensions under consideration 
[1]. Coordinate systems other than Cartesian, more than one spatial 
dimension, and the boundary conditions all can pose problems that are 
either extremely difficult or Impossible, to solve by analytical 
methods. Numerical methods thusly become the only possible solution 
method If the problem complexity Is not to be compromised. Typically 
the ability of a particular method to predict a field variable Is 
tested by numerically solving a problem for which a known exact 
solution Is available. The ability of the method to predict the exact 
results Is a measure of the confidence that can placed In a solution 
where no exact solution exists or experimental test results are 
unavailable. 

The objective of the work to be presented Is to extend, expand, 
and compare an explicit exponential finite difference technique first 
proposed by Bhattacharya [2]. To date the method has only been used 
for one-dimensional unsteady-state, heat transfer problems In 
Cartesian coordinates. 

The method has been expanded, In this report, to allow application 
to a variety of problems. The exponential method will be extended here 
to the case of one-dimensional unsteady heat transfer In cylindrical 
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coordinates. Also, It was used In two- and three-dimensional unsteady 
heat condition. Other cases of Interest that were solved using this 
approach Include temperature varying conductivity In one-dimensional 
heat transfer and the development of the temperature field In laminar 
Couette flow. Solutions of the above cases were either compared to 
exact solutions or to results obtained by alternative numerical 
techniques . 

One final application of the exponential finite difference 
algorithm was made for nonlinear partial differential equations. 
Burger's equation along with the boundary layer equations are solved 
using the exponential method. Thus, a demonstration of how to apply 
the method to nonlinear problems Is described. 

The results of all the different cases considered In this study 
Indicated that the exponential technique produced results that were 
more accurate than those found through other numerical techniques. All 
exponential computer codes and those of the other competing numerical 
analysis were run In double precision on either the IBM-3033 or the 
Cray X-MP mainframes. The computer codes developed for the exponential 
finite difference method and other numerical techniques used for 
comparison are contained In the appendix of this report. 


II. - ANALYSIS 


The Exponential Finite Difference Algorithm 


An explicit exponential finite difference algorithm as first 
derived by Bhattacharya [2] will now be presented. The method can be 
applied to many of the partial differential equations found In 
engineering and physics. The diffusion equation as It applies to 
conduction heat transfer will be used In the demonstration that 
follows. In reference [2-3] the method was derived for one-dimensional 
conduction heat transfer In Cartesian coordinates. To show how the 
method can be extended, a derivation parallel to the one presented In 
reference [2] will be made for unsteady state heat conduction In 
two-dimensions. Equations of this type are typically solved 
numerically by a variety of methods [4]. 

For two-dimensional heat transfer In Cartesian coordinates with 
constant material properties, the appropriate partial differential 
equation Is [5]: 


To Initiate the exponential method a product solution Is assumed and Is 
written as: 



(1) 


T(x,y,t) = 4 >(x)n»(y)e(t) 


( 2 ) 


The Initial conditions of the problem are assumed to be 


T(x,y,0) = f ( x , y ) 
e(0) = 1 


(3) 
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Now taking the appropriate devlatlves of Eq. (2) with respect to the 
Independent variables .yields 


3T 


36 


3t = dt ; 


3 2 T 


= i ^6 


. 

.2 ’ 


3 2 T 


= <$>e 


3x 3x" 3y 

Substituting Eq. (4) Into Eq. (1) produces: 

£ 

3x~ dy 

Dividing both sides of the above by gives: 


a 2 * 

ay 2 


(4) 


= <» 1 4 ' e ^~2 + 4 >© ^~2 


1 ae 

e at 


(l3J! + l a!tl. 

I* 3x 2 ^ 3y 2 J " 


(5) 


It can be seen that the variables have been separated. Consequently 
both sides of Eq. (5) must equal a constant, say, k. 

Now examining only the left hand side of Eq. (5), 

1 36 

6 3t " " K 


Multiplying the left-hand side of this equation by 4^/44 gives: 

4>vU 30 

6 <M dt = ‘ K 

which can be rewritten from Eq. (4) as: 

1 3T 

T(x,y,t) 3t = " K 
Direct Integration produces: 

T = c 2 exp j- k 1 1 


Next, the Initial consltlon Is used to evaluate the Integration constant 
giving c ? = T(x,y,0); thus, 

T(x,y,t) = T(x,y,0) exp j- Ktj ( 6 ) 

Returning to Eq. (5) only this time, the right-hand side Is examined; 
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1 , 1 I 

* ax 2 ^ ay 2 j 


Multiply this by and obtain 

v<p\(f 


I ®«* *x 2 


2 2 

eu> a <t> x ea> a or 

ay 2 


or 


a 

0<t>4/ 


64» ^ + e<t> 


Ia 

ax^ ay' 

Equations (2) and (4) are now used which results In: 


a ( a£i 3 T ) 

Max 2 ay 2 j = 


(7) 


The temperature appearing as a coefficient Is replaced using Its 
Initial value so that, 

. 2 , 


T(x,y,0) 


+ a_T i 

ax 2 ay 2 ) 


(8) 


The partial derivative terms can be written In central difference form 
about a node (1,j) as [6]: 

2 T n + T n - 2T n 

afl „ 1 + l.-l 1-1.3 JAA 

2 ~ 2 
ax ( ax) ^ 


( 9 ) 


2 T° + T° - 2T° 
81 „ 'l.Hl 1 1.1-1 *'l.,1 


ay 


Thus Eq. (8) becomes: 


T n + T n - 2T n 

hi j i-i, j jaa 

(ax ) 2 


u 


(ay)‘ 


T n T n 

1.H1 * 'l,J-l 


2T 


Ul = - 


(ay)' 


( 10 ) 
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Now Eq. (10) Is substituted to replace the constant k In the 
exponential of Eq. (6). Making the appropriate substitutions results 
In: 


T 1 J = T ?J 6XP 


aAt 


'l.J L 


T n + T n - 2T n 

l+l, J 1-1 J tllA 

(AX ) 2 


T n T n 

1.1+1 " 1.J-1 


2T 


ixl 


(ay) 

If the grid spacing Is constant (Ax = Ay), then the above equation 
can be written as: 


T n+1 - T n exp |-^t 
i.J i.J ■ ( Ax j 


•y-n , . T n T n * 

T i+u T i-u T u+i T u-i ~ 4T u 


i.j 


(11) 


Note the At that appears In the above Is the time elapsed between 
time steps n and n+1. The temperatures on the right-hand side are 

i. L_ 

the four nearest neighbors to the 1,j node (see Fig. 1). 

In keeping with the notation derived by Bhattacharya [2], the term 

aAt 


fi = 


(ax)' 


( 12 ) 


Is called the dimensionless time step and 


I.J 


,n T n T n T n . T n 

1+1, J * T 1-1,J * T 1,.1+l * T U-1 ~ 4T 1,J 

T n 

I.J 


(13) 


Is called the dimensionless drive number. Thus, Eq. (11) can be 
rewritten rather compactly as: 

T u - T ".J “pH.j} (14) 

It was found [2] that an Improvement In the solution at the n+1 
time step at node (1,J) can be made If the time step Is divided Into 
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a number of equal length time sub-intervals. The method of time step 
sub-intervals can be described as follows. Let us assume that the time 
Interval will be divided Into three Intervals Including the one at the 
end of the time step (Fig. 2). Now returning to Eq. (14), and 
evaluating at the n+1/3 time step, results In: 


T n+l/3 T n 

T i.j - T i.j exp 




T n+2/3 T n+l/3 

T i,j ' T u exp 


Proceeding from the n+1/3 time step to the n+2/3 time step: 

/ Q M n+l/3 \ 

(3 M I,J / 

Finally, the temperature can be written at the n+1 time step: 

T n+1 T n+2/3 M n+2/3\ 

eXP \3 N 1,J ) 


(15) 


(16) 


(17) 


Now substituting Eq. (15) Into Eq. (16) and substitution of Eq. (16) 
Into Eq. (17), produces: 


T n+1 T f « u n ) IQ u n+l/3l f S2 u n+2/3\ 

T U 3 T 1.J eXP (j M t.j| exp {3 H 1.j ) exp { 3 M 1 , J ) 


(18) 


or 



where the M's are then evaluated at the sub-time Intervals and then 
summed for calculation of "T" at the n+1 time step. Equation (19) 
can be written more generally as: 




exp 


m 

ft \ s M n+p/(m+l) 
m + 1 Z-f 
P=0 


( 20 ) 


In reference [2], It was shown that for heat transfer applications, the 
time step can be subdivided as follows: 
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m 



finite heat transfer coefficient 
with flanking nodes outside 
calculation domain. 


heat transfer coefficient -» ® 


(21) 


N = number of nodes In one of the coordinate directions. 

The procedure necessary to determine the dimensionless drive 

numbers will now be described. Since the method Is an explicit 

technique, all information Is known from the previous time step or from 

the previous time sub-interval. The effort Is the centered around 

calculation of the drive numbers at the requested number of time 

sub-intervals for each of the spatial positions (nodes). 

The calculation procedure Is shown In Fig. (3). Because the drive 

numbers are evaluated at sub-time Intervals, the temperature (or any 

other field variable) must also be known at these sub-time Intervals. 

Therefore, the temperature field Is calculated at each sub-time 

Interval, and In turn Is used to calculate the next drive number. The 

% 

drive numbers for each node are summed for all the sub-time Interval 
steps and then used to predict the field variable at the next complete 
time step. This results In a computer storage requirement of 4(N), 
where N Is the number of nodes. 


In the previous section, the exponential finite difference 
technique was extended to two-dimensions. Now the procedure will be 
considered In another coordinate system. In particular the method will 
next be applied to radial one-dimensional, unsteady heat transfer. 

The governing unsteady diffusion equation for a material with 


Extension To One-Dimensional Cylindrical Coordinates 


constant properties Is: 
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*p S - * [* fc (' £)] 

f a 2 T 1 ail 

a L^ + F57 J 


( 22 ) 


31 

at 


Assume that the temperature can be written as the product: 

T(r,t) = ♦(r)e(t) 

The Initial conditions are 

T( r ,0) = f ( r ) ; 0(0) = 1 

Now taking the appropriate derivatives of Eq. (23); results In 


3J 

3t 


<*> 


2 2 

36 aj _ a$ ajr _ a_$ 

at • ar ' 8 ar ’ „ 2 3 ® 2 

ar ar 


(23) 


(24) 


(25) 


Substituting Eq. (25) Into Eq. (23) yields: 

.2 


. 30 

* at = “ 


e 3_$ 1 e M 

ar 2 r 


Dividing both sides by 4 >e brings: 


1 3i 
e at 


i_ a <£ + l a<t> 

4> ar 2 r<|> ar 


= - K 


(26) 


It can be seen that the variables have been separated. Thus, both 
sides of Eq. (26) must equal a constant, k . Now using the time 
dependent side of the above , l.e., 

1 ae 


e at 


= - K 


and multiplying this by © 4 >/© 4 > produces: 

1 31 

T at = K 


Integrating for a particular value of the radial position "r" gives: 

T = C-) exp {- *t} 
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Next the Initial condition Is used and the equation can be written as: 

T(r,t) = T( r ,0) exp {- «t} (27) 

Returning to Eq. (26) and using the radial side of the equation, we have 

2 

1 + 1_ 

4> . 2 r<j> 3 r 
3r , 

Multiplying by e<t>/e<t>, the equation can be rewritten as: 


1 a T 1 3T 

T ar 2 Tr 3r 


Incorporating the Initial condition, 


(28) 


a 3 2 T 1 all 

T(r,0)[ ar 2 r arj 


(29) 


Next the partial derivatives are replaced using central finite 
differences [4]: 


2 j n + j n 2T n 

Vl 1-1 1 


3r 


r n 


(Ar)' 
r n 


3T T H1 T 1-l 


3r 2 Ar 

Substituting Eq. (30) Into Eq. (29): 

T ni * T ?-! - i_ 

1L V (»r) 2 / 


n n 

'HI 1-1 

2 Ar 


(30) 


(31) 


Equation (31) Is used to replace the constant, - k In Eq. (27), l.e., 


r n+l 


t" exp 


At 


T n T n 
T H1 + T 1-l 


- 2T 


(Ar)‘ 


' n _ n 

'hi Vl 

2 Ar 
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Rearranging this brings: 


r n+l 


T” exp 


a At 
( Ar ) 2 


1+1 


+ T 


1-1 



( 32 ) 


Equation (32) states that the temperature at the 1^ radial 
position at the (n+1) time step Is found from the product of 
temperature at 1 th position, n th time step and the exponential 
term that Is composed of a dimensionless time step: 




a At 

Ur) 2 


and a dimensionless drive number: 




T n T n 

T 1+l * T 1-l 


rO 




(33) 


(34) 


It should be noted that this drive number applies to the Interior nodes 
(r^ * 0). For the node at r = 0 the dimensionless drive number 
becomes: 


M 


n 

1 



(35) 


Finally Eq. (32) can be written In a compact form as 


T n+1 T n 
T 1 = T 1 exp 


HI 


(36) 


The sub-time Interval evaluation for Eq. (36) Is the same as that found 
In the two-dimensional Cartesian form as shown earlier. 


Stability of the Exponential Finite Difference Method 
With few exceptions, explicit finite difference procedures for 
solving partial differential equations are Inherently unstable, unless 
certain numerical condltons are satisfied. These conditions take the 
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form of a grid size and/or time step requirement written In terms of 
parameters of the given problem. If these stability conditions are not 
met, the solution will diverge, often rather drastically. On the other 
hand, the stability requirements can make explicit methods Impractical 
for a particular application by requiring an unrealistically small grid 
or time step. Nonetheless, these conditions must be known prior to the 
use of any explicit finite difference procedure. 

There are a variety of methods that have been used to establish 
the stability constraints of a finite difference procedure: some are 

very elementary, some quite Involved. In essence, the methods seek to 
find an expression for the amplification factor which Is the ratio of 
the current solution result to that In the previous step. If the 
absolute value of the ratio Is less than one the method Is stable. 
Determination of the amplification factor for the exponential finite 
difference method Is particularly convenient, as has been shown In 
[2]. For the one-dimensional cylindrical coordinate case, the 
amplification factor E can be readily defined as 



(37) 


or from Eq. (32) for an Interior node: 

T n+1 


l = 


±_ 


exp 


K} 


( 38 ) 


So for stability to exist as at and ar approach zero: 
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11m 1 1 1 < 1 

At -» 0 
Ar -» 0 


(39) 


To satisfy this requirement the exponent In the eponentlal of 
Eq. (38) must obviously be less than or equal to zero. Since the 
components that make up Q In that exponent are all positive, this 
Implies that the dimensionless drive number will dictate whether or not 
the stability criteria Is met. For the cylindrical coordinate case the 
dimensionless drive number must satisfy: 

„ T n\ 


H n | f T Ul * T 1-l 
1 T ; 


Multiplying by this becomes: 



- C , 
* 


< 0 


(40) 


( T n + T n _ 2T n \ + AE_ / T n - T n ) < 0 

\ 1+1 1-1 1 J 2 r 1 ^ 1+1 1 1-1 J - U 


Define 13 = Ar/2r^ and rearrange to get 


(1 + 0 ) T " +1 + (1 


0)T"_-, < 2T" 


t" > \ (1 + 0)t" +1 + (1 - 


(41) 


Equation (41) needs to be satisfied otherwise an unstable 
condition can exist. In Eq. (41) as Ar -» 0, or equivalently 0 -* 0, 
the stability condition becomes: 

..n 


T n ^ I T n T , 

T 1 - 2 T 1+l + T 1-l 


(42) 


Also for one-dimensional cylindrical coordinate case, the node at 
r = 0 has a different stability criteria because the node at 1+1 
does not exist. Since the radial derivative of the field variable must 


equal zero at r = 0. In finite difference terms this can be realized 
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by requiring that t" + 1 = T I J_ 1 at this node. Thus, at the 
origin, Eq. (42) becomes 

T 1 £ T ?_i («) 

As stated by Bhattacharya [2], the dimensionless drive number Is 
the determining factor whether or not the stability criteria Is met. 
However, the dimensionless time. If made large enough, could cause the 
solution to become unstable. Since time sub-interval division Is used, 
the total dimensionless time step Q could become quite large. From 
[2] It was recommended that the dimensionless time step be chosen to 
satisfy the following: 


where m Is the number of time step sub-intervals Involved In the 
calculations. If Q = 5, for example, then m would have to be 
greater than or equal to 9 l.e., nine sub-intervals would be required. 
From Eq. (21) with Infinite heat transfer coefflcent, a minimum of 20 
nodes would be needed. 

Another useful comparison to be made Is the one-node model as used 
In Refs. [2] and [7] where the value of the dependent variable at the 
surrounding nodes Is set equal to zero. For the one-node model, as 
stated In Ref. [7], the exponential finite difference and exact are the 
same (Fig. 4). This figure Indicates that the exponential solution 


remains stable as Q Is Increased. 


15 


Effect of Initial and Boundary Conditions on the Exponential 
Finite Difference Method 

In this section the effect of boundary conditions on the 
exponential finite difference technique will be Investigated. Boundary 
conditions that are typical of heat transfer applications will be 
considered. The conditions to be presented are: (1) finite heat 

transfer coefficient (mixed condition), (2) Infinite heat transfer 
coefficient (Dlrlchlet condition), (3) constant heat flux (Neumann 
condition), and (4) time varying. Initial conditions, where the field 
variable Is equal to zero, will also be discussed. 

Finite Heat Transfer Coefficient 

For this boundary condition, the method used In reference [2] will 
be utilized. Numerical Implementation of the boundary condition 
requires that a node be placed outside the solid In the surrounding 
medium. This external node will be used In the finite difference 
equation at the solid surface. One-dimensional heat transfer In a 
cylinder (T = T(r,t)) will be used to demonstrate the procedure. 

Using the exponential finite difference method for T a (r,t), it 
was found earlier (Eq. (36)) that 

T 1 = T 1 exp |QM 1 J 


(45) 
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where the dimensionless time Is given by 

a At 

Q = 2 

(Ar)' 

and the dimensionless drive number by: 


j n + y n _ 2T n \ / T - T 

1+1 1-1 1 \ . ar / 1-1 1+1 


'1 




2r 


1 




( 46 ) 


The thermal condition at the surface Is found by equating the 
conductive and convective heat fluxes at the surface: 




- h ( T | 

1 » ' T ") 

r=R 

V 1 

1 r =R / 


(47) 


Using the node numbering as shown In Fig. 5 and a central difference to 
express the deveratlve, Eq. (47) becomes: 


(y ) 


Solving for Tq, the temperature of the external node, yields 


,.n 2h ar 
T 2 * — 




(48) 


(49) 


Let B = h ar/k (Blot number [8]), thus Eq. (49) Is written: 

T n 0 - t!| + 2BT" - 2BT^ (50) 

Now that an expression for the temperature at the external node Is 
known. It can be used In the expression for the temperature at the 
surface node. This results In: 


T n+1 x n / _ u n 1 
T^ = T-j exp js H| j 


where 


j 

T 2 - 2T " * 2BT . * T 2 - 2BT " 

. Ar 

\^: - t?) 

= [ 

J 

+ R 



(51) 
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and the drive number can be further simplified to: 


1 

2 Tj - (1 + B)T" + BT ro 

ArB 

[t. - 

J 

( 

J 

+ T" 

T ? J 

j 


Infinite heat transfer coefficient 
If the heat transfer coefficient In Eq. (47) Is placed on the 
left-hand side of the equation and then allowed to become Infinite, It 
Is seen that the surface temperature will equal the temperature of the 
surroundings. Thus, the boundary condition In which h -» <® Is 
Identical to that In which a boundary temperature Is held constant. In 
the calculation procedure, these Isothermal boundary nodes are only 
needed for calculation of the temperature field at the surrounding 
nodes. For example In a two-dimensional square grid with a total of 
121 nodes calculation would be reduced to a total of 81 nodes If the ; 
temperature Is specified for all four boundaries. 


Constant Heat Flux 

For a constant heat flux applied to the boundary surface, the same 
procedure as was used In the finite heat transfer coefficient case will 
be utilized. The condition at the surface for one-dimensional 
cylindrical coordinates Is given by: 

(53) 

r=R 



An external node Is placed outside the solid as was done earlier. 
Equation (53) can then be written 



q = k 


(54) 
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Solving for the external node temperature results In: 


T n 2q Ar T n 
'O " k '2 


( 55 ) 


At the surface the exponential finite difference equation Is: 


T n+1 T n 

= T-j exp 


ry n + j n _ 2T n \ /T n - T n 
*0 '2 _ 1 ) ar f !o 12 

T n / + 2R V ,n 


1 


1 


and substituting in Eq. (55): 

2q Ar T n 
k '2 


Rearranging this produces: 


t" +1 = t" exp 


n 




at 

2R 



+ T, 




= exp 


Q 


2 


q Ar 
k 




+ 


(Ar) 2 q 
Rkf" 


(56) 


Time Varlng Boundary Conditions 

This condition Is similar to the constant boundary temperature 
condition except that the boundary temperature must be Incremented as 
the calculation marches In time. The boundary condition must reside In 
the time step loop which Is shown as the outer most loop In Fig. 3. 

The temperatures on these boundaries are Incremented and held constant 
as the subtime Interval calculations are made. 

Dependent Variable Initially Equal to Zero 
One last condition that can exist wherever there Is Initial zero 
temperature (T(x,0) = 0) needs to be discussed. If this condition Is 
encountered, then the following substitution should be made or else the 
exponential finite difference method will not work. This can be 
readily seen by examlng any of the numerical equations e.g., Eq. (56). 
Since the Initial temperature would appear In the denominator of the 
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exponent In the exponential, problems would ensue. To circumvent this 
difficulty define a new variable T, such that T(x,t) = 1.0 - T(x,t). 
Now the exponential finite difference equations described above can be 
utilized by simple replacement of the T variable with the T 
variable. 


III. NUMERICAL COMPARISON OF THE EXPONENTIAL FINITE DIFFERENCE 


METHOD TO EXACT SOLUTIONS AND OTHER NUMERICAL TECHNIQUES 
The final product of any numerical study Is how well the given 
method performs when compared to known exact solutions or to other 
numerical techniques. The exponential finite difference method will 
now be applied to the following cases to demonstrate the capability of 
the method to solve the diffusion equation: 

(I) One-dimensional heat conduction In cylindrical coordinates with an 
Infinite and a finite heat transfer coefficient at the surface, 
unsteady state, 

(II) Two-dimensional unsteady state heat conduction In Cartesian 
coordinates , 

(III) Solution of Laplaces equation, Cartesian coordinates, 

(1v) One-dimensional heat conduction In Cartesian coordinates, unsteady 
state with temperature varying thermal conductivity, 

(v) Steady state Couette flow, 

(vl) Three-dimensional heat conduction In Cartesian coordinates, 
unsteady state. 

One-Dimensional Heat Conduction In Cylindrical Coordinates 
The one-dimensional cylindrical coordinate heat conduction case 
with temperature as a function of time and radial position will be 
Investigated for Infinite and finite heat transfer coefficient. The 
exact results for both cases can be found In Ref. [9], 
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For Infinite heat transfer coefficient on the boundary surface the 
exact result Is given In [9] as: 


T(r,t) 


T n - T 
U 00 


= 2 


-axjlt 

_ e w> 

m=l (V) V X m R) 


E 


(57) 


where X-R Is the m th zero of 
m 

J 0 (X-R) =0 (58) 

The results of both the exact analysis and the exponential finite 
difference method are shown In Table I. As can be seen from the 
tabulated results, exponential finite difference results approach the 
exact solution as the number of nodes Is Increased or as the 
dimensionless time step Is decreased. 

When the heat transfer coefficient has a finite value at the 
surface, the exact solution from [9] Is: 


T (r,t) - T a 
T n - T 

(J 00 


= 2B 


E 


m=l 


-ax|t 

e J 0 (X-r) 


2„2 


(x? R 
\ m 


+ B‘ 


'>W> 


(59) 


where B = hR/k (Blot number) and X- (characteristics values) are 
given by (for cooling) : 

<V»W> - B W> ■ 0 (60) 

The results are shown In Table II for various values of the Blot 
number. As would be expected the solution approaches the exact 
solution as the number of nodes In Increased. The size of the Blot 
number did not seem to effect the accuracy of the solution. As the 
elapsed time of the solution proceeded, temperatures predicted by the 
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exponential finite difference method approached the exact result. Also 
the results Indicated that reducing the size of the time sub-interval 
Increased the method's accuracy. 

One last comparison will be made while Investigating the 
exponential finite difference technique In one-dimensional cylindrical 
coordinates. The problem situation Is shown In Fig. (6) and applies to 
a cylindrical annulus with the following Initial and boundary 
conditions : 

T(r.O) = 0 

T(R 2 .t) = 1.0 (61) 

fr < R 1 ’*> 

In Ref. [4] this problem was solved numerically using a 
characteristic-value solution. A comparison of results Is shown In 
Table III for the exponential method using the same grid spacing as In 
[4] and for the case where grid spacing Is halved. The results are 
seen to compare quite well with the finer mesh being slightly closer to 
the value from Ref. [4] especially during the first few time steps of 
the solution. 

Two-Dimensional Heat Conduction In Cartesian Coordinates 

The exponential finite difference technique will now be applied to 
the two-dimensional heat conduction problem In Cartesian coordinates 
shown In Fig. (7). The exponential finite difference method will be 
compared to the solution of this problem, as performed In Ref. [4], 
using the alternating direction Implicit technique (ADI). The results 
of the two numerical techniques and the exact analysis are shown In 
Table IV. The temperature Indicated for comparison Is that at the 
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origin x = y = 0, shown In Fig. (7). As maybe seen, the ADI technique 
does not predict the temperature as accurately as the exponential 
finite difference method at the first time value shown In Table IV. 
However, as the amount of elapsed time Increases either method does a 
very good job at predicting this temperature. When the number of grid 
points was Increased, by halving the spatial Intervals, the exponential 
finite difference method was found to be more accurate for all the time 
steps . 

Since the ADI method Is one that requires simultaneous solution of 
equations In the two coordinate directions, the time step size can be 
made large. The exponential finite difference technique must have the 
dimensionless time step kept below 0.25 to keep the solution stable. 

So the required CPU time for the exponential method Is higher for this 
application. 


Solution of Laplaces Equation 

Since the exponential finite difference method has been used for 
two-dimensional unsteady state conduction, a natural extension with 
little additional effort would be to use this method to solve Laplace's 
equation. This can be Implemented In the exponential finite difference 
method by just allowing the solution to march In time until no further 
change In the field variable Is Indicated. 

As an example, the problem as shown In Fig. (8) will be solved and 
the results compared to those given In Ref. [4], In the referenced 
work, the solution was found by using a Gauss-Seldel Iterative 
technique. 
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A comparison of results along a diagonal from the position (x = 0, 
y = 1) to (x = 1, y = 0) Is presented In Table V for two different grid 
spaclngs. As can be seen, the solutions are nearly Identical with the 
exponential method requiring a smaller number of Iterations (or time 
step Increments) to reach a slmlllar result. 

One-Dimensional Unsteady State Conduction With Temperature 


The effect of temperature varlng thermal conductivity will now be 
Investigated using three different numerical schemes: a pure explicit, 

the exponential method and an Implicit technique. The problem to be 
solved Is Illustrated In Fig. (9a). The thermal conductivity as shown 
In Fig. (9b) Is assumed to be a linear function of temperature. 

The exponential finite difference method will be applied first to 
the given problem. The governing partial differential equation Is [1]: 


From [1], Eq. (62) can be changed to a simpler form by using a new 
variable e (the Klrchoff transformation) given by: 


Varlng Thermal Conductivity 



(62) 



R 


k(T)dT 


(63) 


where k D Is the conductivity at temeprature T , and 



il _ ae 
at ” k at 


( 64 ) 


ae = k_ aT 
ax “ k R ax 


Substituting Eq. (64) Into (62) gives: 
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pC p k R /ae\ _ a_ ( v ae\ 
k vat/ ax ^ R ax J 


p c r 


(ft)- 


a 2 © 


ax 


(65) 


Since It has been assumed that the thermal conductivity Is a linear 
function of temperature, 


k(T) = k R (l + I3T) 


( 66 ) 


Now returning to Eq. (63) and substituting In Eq. (66), we have: 


6 



+ (3k r T ) dT 


Direct Integration yields: 

e = (T - T R ) {l + § (T + T r )} (67) 

Equation (67) provides the relationship between the variable T and 
the new variable e. 

Returning to Eq. (65), 
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at 


k 

pC 


a 2 e 
p ax 2 


( 68 ) 


Equation (68) Is In a form now that the exponential finite 
difference can be applied. The resulting equation In the variable can 
be shown to be given by: 


e^ 1 = e^ 1 exp 


At 


P C p (Ax)‘ 




+ 0 


1-1 
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?) 


(69) 
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Evaluating Eq. (67) at node 1 brings 


or 


- ft - t r) I’ * § ft * t r) 
9 ? ’ ft - t r) * ! (ft) 2 * t r) 


( 70 ) 


Substitution of Eq. (70) Into Eq. (69) at the appropriate time steps 
and nodal locations will give: 


ft 


exp 


n+1 





ft - t r) * 1 

ft) 2 ft 

’(ftl * ftl - ft ) * II 

[(ft,f * ftj 

! - *(ft 2 l 

ft - ft * ! [(ft 2 - 1] 


(71) 


where 


Y = 


at 


pC p (ax) 

If Tr = To, = 0.0, Eq. (71) becomes: 

T r * ! ft *') 2 ■ ft * i ft ■ 

ft.,* ft,)-* T? - ir^.n) 2 - K.,) 2 - Z (T?) ; 


exp 




* |(T?) 2 


(72) 


r n+l 


The equation for Is a quadratic with the right-hand side of the 

equation all being known at time step n, so define a variable 
such that 

ft.l * ft) - 2T ") t |[( T U!) 2 * (ft))' - ? W) 2 ] , 

* I ft) 2 


■ 

I 
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Equation (72) then becomes: 



( 74 ) 


Solving this and using the positive root results In: 


n+1 _ 1 
1 = B 


- 1 + yi + 2^0 


(75) 


where 


0 > 0 


Equation (75) In conjunction with Eq. (74) are Implemented In the 
exponential finite difference solution sequence. In this case the 
conductivity as well as the temperature field must be kept track of on 
the sub-time Interval level. The dimensionless time step, Q, and the 
rate of conductivity change, 0, must be both considered when choosing 
the step size so the solution does not become unstable. For this 
method, the term (yk^/m+l) in the exponential was considered at 
Its maximum possible value and the time step was adjusted to retain 
stability. This criteria was chosen so that 


The next method to be Investigated for the temperature varying 
conductivity problem will be the pure explicit method. As stated 
earlier the governing partial differential equation for one-dimensional 
conduction Is given by: 


Using the chain rule this equation can be put Into a nonconservative 



(76) 


form. 
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r aT . a 2 T A ak /aT x2 

it - * ^ ♦ ST (ix 


( 77 ) 


Assuming the same linear profile as before: 


k(T) = k R + 6k R T 


then 


(78) 


— - 8k 
aT ,3K R 

Substituting Eq. (78) Into (77) will give: 

0 C P H ■ k R<> * BT) ^ * 0k R («)* 


or 


ai 1 1 nT , 3?i n /ai\ 
3t - a R (l + 8T) 2 ♦ 0<* R ^- X j 


(78) 


where 


a R " P C r 


Using central space differences and a forward time difference we may 
write: 


j n - j n 
aT HI 'l-l 


ax 


2AX 


T n+1 T n 
aT 'l " 'l 


at 


At 

r n 


2 T n + j n 2T n 

all HI 1-1 1 

? ~ ? 
ax (ax) 

Substituting Eq. (79) Into (78) produces: 


n+1 T n 

1 - pF- 1 - «*(* * bt ?) 


T n T n OT n 
T U1 * T 1-1 - 2T 1 

(AX) 2 


( 79 ) 


+ Ba r 


1+1 


- T 
2bx 


1-1 
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This may be simplified and written as: 



( 80 ) 


where 



Equation (80) can be used to directly solve for the temperature 
field at the next time step. Some additional care must be used to keep 
the solution stable as the size of the dimensionless time step Q, 
and the rate of conductivity change, (3 both effect the solution. 

The final method to be Implemented for comparative purposes Is the 
Implicit method. Starting with Eq. (78), we have: 


To avoid a solution sequence that would require the solution of 
nonlinear algerbralc equations, the following will be assumed; 

(a) The term (1 + 0T) can be replaced with (1 + BT 1 ^) 

(b) The squared first derivative term can be replaced by 



This Is a linearizing technique known as lagging the coefficients. 
Substituting these Into Eq. (81) will produce: 



(81) 



Further simplification gives: 
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T? +1 - T" = ft (l + 6T 


1 


where 


Now define: 


?)( T ? 


n+1 + T n+1 - 2T 
+1 1 1 -1 


r) 


flB /,n T n \ /n+1 T n+1\ 

4 \'i+i " 1 i-i j \'m ~ 'i-v 


Q = 


At 


(AX)‘ 


n (3 / T n 


XV = 


- /j n - T n \ 
4 [ HI 1-1 j 


II i OT ll 

n^ = 1 + 0T^ 


Substituting Eq. (83) Into (82) yields: 


( 82 ) 


(83) 


r n+l 

'l 


T n / n ) T n+1 T n+1 0T n+l 
- T 1 - a ("I T U1 * T 1-l - 2T 1 


* t"* 1 

*i ri+i 


r n+l 
1 1-1 


Simplifying this results In: 

T n+1 
'l 


L ^ n \ T n+ 1 n „ n\ T n+1 n rt n\ T n 

y 2ftn i j T i+i (^ Qn i + flX i J ~ T i_i (^i - ftx i J = T i 


(84) 


The equation shown above Is now In a form that can be used In the 
Thomas Algorithm [4], l.e., Eq. (84), can be written as: 


where 


T n+1 . _n+l T n+1 _n 

a 1 T 1_l + b i T i + c i T i+i = T i 


a 1 = “ Q ( n 1 " X l) 


b 1 = 1 + 2Qn ^ 


C 1 = 




/ n ^ n \ 

("i * N ) 


(85) 


Equation (85) can now be solved using a tri-diagonal matrix routine. 
The variables a^, b^ and c^ must be evaluated at each 
position and time step as their values change as the field variable T 
changes. 
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A comparison of results of the three methods can be found In 
Fig. (10) and Table VI. Figure (10) shows the temperature field 
through the slab cross-section. From this. It Is evident that the 
exponential and pure explicit methods give very slmlllar results. The 
Implicit method predicted higher temperatures closer to the slab 
surface and lower temperature at the slab centerline then either of the 
two explicit methods. In Table VI the results at the slab center are 
shown for various elapsed times. As can be seen, all three methods 
agreed with each other to within a few percent. 


The Steady State Temperature Field of a Couette Flow 
Another application of the exponential finite difference method 
will now be presented. The problem to be Investigated Is the 
developing temperature field In laminar Couette flow [7]. The problem 
statement Is Illustrated In Fig. (11). Neglecting viscous effects, the 
governing equation Is given by [10]: 


Neglecting conduction In the x-dlrectlon or assuming that the 
convection term Is much greater than the conduction term, Eq. (86) 
becomes : 


U 


al 
x ax 


a 



From Fig. (11) using the expression for the velocity In the 
x-dlrectlon, U x = Uy/L, Eq. (87) becomes 


( 87 ) 
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( 88 ) 


3x " Uy ay 2 


Equation (88) Is now In a form where separation of variables can 
be Implemented. Following the same procedure as Indicated earlier to 
find the exponential finite difference equation. It can be shown that: 


The procedure utilized here Is that the solution marches In the 
x-dlrectlon Instead of time as was the case for the previous examples. 
Information from the last x-posltlon and y-dlrectlon are used to 
determine the dependent variable at the next x-posltlon. 

Results of Implementing this method are shown In Fig. (12). The 
temperature field Is shown for three x-locatlons for two different 
values of L/U. The results Indicate that as the upper plate velocity 
U Is Increased, the propogatlon of the temperature change In the 
y-dlrectlon Is slowed down. 

Unsteady State Heat Conduction In Three-Dimensional Coordinates 

The final application of the exponential finite difference method 
to the diffusion equation will be that of three-dimensional, unsteady 
state heat conduction. The exponential method, a pure explicit method, 
and an Implicit method (method of Douglas, [11]) will be compared to 
exact solution for the situation shown In Fig. (13). 

The exact solution to the problem Illustrated In Fig. 13 Is given 



(89) 


In Ref. [10] as: 
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00 00 00 



where a, b, and c are the widths of the cube In the x, y, and z 
directions respectively. Equation (90) will be used to determine how 
well the numerical techniques predict the temperature distribution. 

The exponential finite difference technique will be Investigated 
first. The sequence to be followed for determining the finite 
difference equation Is the same as presented for the earlier cases. 

The procedure for this three-dimensional case consists of the following 
stepped procedure: 

(1) Linearize the partial differential equation 

(2) Assume a product solution 

(3) Separate time from spatial dependence 

(4) Solve for time dependence 

(5) Insert the appropriate spatial finite differences Into 
exponential term that results from step 3 

Based on this procedure the three-dimensional exponential finite 
difference equation can be shown to be: 
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becomes : 


T 


n+1 

U,k 



exp 


m 

ft \ N M n+p/(m+l) 
m + 1 1,J,k 

P=0 


(92) 


where m Is the number of subtime Intervals, ft Is the 

dimensionless time step, and m" j k the dimensionless drive number 
given by. 



Equation (92) will be used for all Interior nodes In Fig. 13. This 
equation, as well as those that result from the other analysis, will be 
adjusted along the Insulated boundaries to take Into account the 
boundary condition that exists there. 

The next method to be applied to this three-dimensional case will 
be the pure explicit method. The finite difference equation for this 
method Is given by [11]: 

T 1,j,k = T 1,j,k (1 " 6Q) Q (Vl.J.k T 1-l,j,k T 1 ,J+1 ,k 

T n T n T n v 

'l.J-l.k + 1 ,J ,k+l 1 ,J ,k-l j 


( 94 ) 
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a. At 

where fl = j ; Ax = Ay = Az 

(A X)' 

As shown In Ref. [11] the dimensionless time step Q must be: 



( 95 ) 


to ensure stability of the method. 

The last numerical technique that will be applied Is the Method of 
Douglas [11]. This method Is Implicit, and the spatial directions are 
considered sequentially In the x, y, and then z directions 
respectively. The Intermediate temperatures U (found from x-dlrect 
sweep) and V (found from y-dlrectlon sweep) are used to calculate the 
actual temperature field variable T (found from z-dlrectlon sweep). 

The equations that are solved sequentially are presented as follows. 



'i.J .fc ~ T ”. J.k . 1 t 2 ( u , T n \.1 4 2 [ v , T n ) 

a At ■ 2 \ p.J.k I.J.kJ 2 y \ I.J.k 1J,k/ 

t 2 / T n 
+ S z ( T U,I 


T n+1 T n 
l.-l.k ~ I.J.k _ 1 ,2 


(97) 
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+ T' 


+ 1 6 2 ( V, 


+ T 


2 x [ IJ.k I.J.k / 2 y \ 1J,k I.J.k 

1 .2 / T n+1 T n 

+ 2 6 z T 1,J,k T 1,J,k 


(98) 


where the finite difference operator In the x-dlrectlon , for example, 
would be: 


6 = 


n T n 

m . J .k + '-1 1 J >k 


2T 




(Ax ) 4 


(99) 


Equations (96), (97), and (98) must be solved successively because the 
variable U Is used In equation (97) to find V and so on. Since the 
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method operates on one spatial direction at a time, the Thomas 
Algorithm can be utilized. In the case of finding the U variable, 
the y and z nodal positions are held constant for all the 

x-dlrectlon nodal positions (Fig. 14). This process Is repeated until 
all y and z nodal values for the x-dlrectlon variable U are 

calculated. This procedure Is then repeated In a similar way for the 
V variable and then finally for the actual temperature field variable. 

The results from the three different, three-dimensional solution 
methods are shown In Table VII. The exponential finite difference 
method described above outperformed the pure explicit and the method of 
Douglas for all positions as shown In Table VII. 

In Ref. [11] nine different methods to solve the diffusion 
equation In th^ee dimensions were Investigated. The method of Douglas 
was the preferred method because of Its accurate results and low 
computer CPU time. In that study the pure explicit method required the 
lowest amount of CPU time with the method of Douglas requiring 
approximately four times as much. In the present study all three 
methods were run on two different mainframe computers to Investigate 
how well these three methods compared In CPU times. The results are 
shown In Table VIII. All three methods were exercised for the same 
number of time steps. As Indicated, the exponential method was 
approximately three times faster than the method of Douglas but still 
slower than pure explicit method. From these results It could be 
concluded that the exponential method would have been chosen as the 
preferred method had It been used In competition with the nine 
numerical methods as described In Ref. [11]. 


IV. - EXTENSION OF THE EXPONENTIAL FINITE DIFFERENCE METHOD 
TO NONLINEAR PARTIAL DIFFERENTIAL EQUATIONS 


The exponential finite difference method will now be applied to 
two different nonlinear problems. The problems to be addressed will be 
the viscous Burger's equation and the boundary layer equations (steady 
state flow over a flat plate). 

The viscous Burger's equation Is given In Ref. [12] as: 


au .. au a 2 u 

at + u di = v 7~2 

ax 


To allow application of the exponential method to 
equation must be first linearized. So letting U 
the nonlinear, term and rearranging the equation, 

au . au a 2 u 


( 100 ) 

this equation, the 
= A = constant, for 
gives : 

(101) 


Assuming a product solution of the form 

U(x,t) = <t>(x)e(t) 

and taking the appropriate derivatives, Eq. (101) 

0 . se . ve J 


becomes 


Division by <|>e gives: 


2 

1 a© A a^ u a <j> 

e at = ” <t> ax + 4 > ax 2 


k = constant 


( 102 ) 


It can be seen that the terms are now separated. As has been shown 
earlier, the left-hand side of Eq. (102) can be written as: 
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U(x,t) = U(x,0) exp {- Kt} (103) 

Now returning to Eq. (102) and examining the x-dependence, we have 


A 3<t> v 3 <t> 

" * dx * ♦ ax 2 ' 

Multiplying both sides by <$>e gives: 




Ae t* + v6 ~ ^ = - k9<}> 

ax ax^ 


This can be written In finite difference form as 
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n I U U1 - U l l 
1 \ 2 AX 
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u 1+l Vl 
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(104) 


This Is used to replace the exponent In Eq. (103), thus 
/ 


u" +1 = u" exp 


Atv 


(AX) . 


- AS. ( u n 

2 v H+l 


n \ A*? i + U 1 

“i-i ) * ( 


i - 2U ? 


(105) 


\ ", 

Equation (105) Is the exponential finite difference equation for the 
viscous Burger's equation. An example will now be used to demonstrate 
the method. 

An exact steady state solution to Burger's equation Is available 
for the following conditions 

U(0,t) = U Q 
U( L,t) = 0 

The steady-state solution was given as [12]: 


U(x) = U Q U 


1 - exp 10 R 



1 |j * exp (o R e (* - 1 )) 
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where 


Re 


L 


V 


(106) 


and U, Is the solution of the equation 


U 1 - 1 ( 

U + -| - exp j - U 1 Re L 


The exponential finite difference method will be now used to 
numerically solve the problem stated above. However, for the stated 
conditions, a problem arises with the portion of the velocity field 
Initially at zero. To overcome this difficulty, the substitution 
method described earlier will be used. A new variable will be defined 
such that 

U = U Q - U 

and Burger's equation then becomes: 


Using Eq. (107) the same method of separation of variables must be 
performed on the U variable. The problem Is now solved for the U 
variable and the substitution shown above Is then made to find the U 
variable. The exponential finite difference equation for U can be 



(107) 


with the following Imposed conditions. If (U Q = 1): 

U(0.t) = 0 


(108) 


U(L.t) = U Q 


shown to be: 
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(109) 


The results obtained by applying Eq. (109) and the conditions In 

Eq. (108) are compared to the steady state exact results of Eq. (106) 

and are shown In Fig. (15). The results from the exponential method 

were nearly the same as the exact method. The exponential method was 

allowed to march In time for quite a number of steps without special 

2 

treatment of the dimensionless group Atv/(Ax) which could have 
been altered to allow convergence to the actual solution In less time 
steps . 

Another application of Burger's equation was made to Investigate 
the effect of the diffusion term. The results for the variation of v 
over four orders of magnitude are shown In Fig. (16) for the same 
Instant In time. At the two lower v values, the total range of the 
field variable takes place over a small number of nodal positions. A 
better approximation could be made for these cases by using a finer 
grid. Also Included on Fig. (16) Is the solution of Burger's equation 
by a pure explicit technique. For the value of » chosen, the 
solution oscillates around the predicted solution found from the 
exponential method. The pure explicit solution was found using the 
same number of nodes and the same dimensionless step size. When the 
solution oscillates, as the pure explicit solution did, the resulting 
velocity field can contain physically Impossible values. 
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The last application to be Investigated will be for the 
development of a laminar boundary layer on a flat plate (Fig. (17)) 
In Ref. [10] the steady state formulation Is given In terms of the 
following three partial differential equations: 
for continuity: 


for momentum: 


for energy: 


3U 3 V n 

« * sy ■ 0 


.. 3U . „ 3U a 2 u 


u ai , v ai , „ a!i 


ax 


ay 


ay 


with the boundary conditions: 

U(x,o) = 0 

V(x,o) = 0 
T(x,o) = 0 


U(o,y) = U ( 

V(x,L) = 0 
T(o,y) = T 


( 110 ) 


( 111 ) 


(112) 


(113) 


v and a are the momentum and thermal dlffuslvltles respectively. 

Equations (111) and (112) can be solved by using the method 
presented for the viscous Burger's equation. The only difference Is 
that the solution will march In the x-dlrectlon Instead of time. 
Keeping this procedure In mind, results of the separation of variables 
for Eqs. (Ill) and (112) were found to be: 
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(114) 


(115) 
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The continuity equation Is written as [12]: 


v 1+ l _ J* 1 

J ' J 


_ay_ 

2 AX 


+ U U1 
j J J-l 


- u 


M 


(116) 


Equations (114) and (115) are first solved using a spatial 
sub-increment as was done for the cases when time was the marching 
direction of the solution. After this, the continuity Eq. (116), Is 
then solved. 

The results of this application are shown In Fig. (18) for a 
Prandtl number equal to 0.72. As can be seen the thermal boundary 
layer was outside the velocity boundary layer as would be expected. 

The results with the Prandtl number equal to 0.72 were compared to the 
exact solution as presented In Ref. [10]. A downstream position was 
chosen and the results are compared In Table IX. The exponential 
method results were In good agreement with the exact results. 


CONCLUDING REMARKS 


In conclusion, an exponential finite difference technique has been 
extended to other coordinates systems and expanded to handle problems 
In two and three dimensions. The method has direct application to 
linear partial differential equations such as the diffusion equation 
and can be extended to solve nonlinear equations. 

The method was applied to the following cases: 

(1) One-dimensional, unsteady state heat transfer In cylindrical 
coordinates, Infinite and finite heat transfer coefficient. 

(2) Two- and three-dimensional, unsteady state heat transfer In 
Cartesian coordinates. 

(3) One dimensional heat transfer, with temperature varying thermal 
conductivity. 

(4) Developing temperature field In laminar Couette flow. 

(5) Nonlinear partial differential equations (Burger's equation 
and boundary layer equations) 

The exponential finite difference method predicted the field 
variable with a higher degree of accuracy In those cases examined where 
the exact solution was available. When extended to three dimensions, 
the. accuracy was still higher for the exponential finite difference but 
the computer CPU time was Increased. When the exponential method was 
compared to other numerical techniques, the results were found to be 
very comparlable. 
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In conclusion, the results predicted for the exponential finite 
difference algorithm for the cases presented In this study demonstrated 
that: 

(1) Field variable was predicted with a higher degree of accuracy 
than other numerical techniques where exact solutions exist. 

(2) The method can be applied to linear and nonlinear partial 
differential equations with dependent variables that can be 
separated . 

(3) The stability of the method Is the same as that of pure 
explicit methods, where the sub-time Interval step size 
determines the stability. 
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TABLE I. - COMPARISON OF RESULTS FOR DIFFERENT DIMENSIONLESS TIME STEP 
FOR ONE-DIMENSIONAL HEAT TRANSFER IN CYLINDRICAL COORDINATES WITH 
INFINITE HEAT TRANSFER COEFFICIENT AT THE SURFACE. INI11AL AND 
BOUNDARY CONDI I IONS ARE: 

[h -» “>, T(r,o) = T .0, T(R,t) =. 0.0, a . a At/(Ar) 2 , o = T.O M 2 /s, 

N = number of nodes, m - number of sub-time Intervals.] 


t. 

sec 

From surface 
r-dlstance 
(M) 

N = 11 
m = 4 
ft = 1 .0 

N * 21 
m = 9 
ft = 1 .0 

N * 21 
m =» 9 
ft * 2.0 

N a 21 
m » 9 
ft a 5.0 

Exact 
analysis 
ref. [9] 

0.1 

0.1 

0.127004 

0.126768 

0.126819 


0.126669 

.1 

1.0 

.062431 

.852204 

.853083 


.848368 

.5 

.1 

.011969 

.011671 

.011680 

.011715 

.011582 

.5 

1.0 

.094334 

.090309 

.090379 

.090652 

.088895 

.5 


Total 

Total 

Total 

Total 




50 steps 

200 steps 

100 steps 

40 steps 




iT^r * °- 2 

= 0.1 

■■ - -- = 0.2 

77 ^- 7 a 0.5 




M f 1 

M + T 

M f 1 

M ♦ 1 



TABLE II. - FINITE HEAT TRANSFER COEFFICIENT 
CYLINDRICAL COORDINATES WITH THE FOLLOWING 


CONDITIONS: T(r,o) » T.O, 

T^ = 0, fi = a At/( Ar ) 2 


Time 

h/k 

R 

Exact 

ref. 

[9] 

Exponential finite 
difference results 

N a 11 
m = 4 
ft a 1 .0 

N a 21 
m a 9 
ft * 5.0 

N a 21 
m » 9 
ft a 1 .0 

0.1 

1 

1 

0.6846 

0.7073 

0.6978 




0 

.9768 

.9814 

.9797 


.2 

1 

1 

.5702 

.5976 

.5857 




0 

.8702 

.8852 

.8780 


.4 

1 

1 

.4132 

.4441 

.4303 




0 

.6420 

.6698 

.6563 


.1 

2 

1 

.5009 

.5285 

.5199 




0 

.9594 

.9670 

.9643 


.1 

5 

1 

.2558 

.2777 


0.2669 



0 

.9265 

.9385 


.9306 
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TABLE III. - COMPARISON OF EXPONEN1 1AL 
FINITE DIFFERENCE METHOD IN ONE- 
DIMENSIONAL CYLINDRICAL COORDINATES 
TO THE RESULTS OF REFERENCE [4]. 

[a = 1.0, At = 1.0 S6C , a At/Ar^ * 1.0. 
N = number of nodes, m = number of 
sub-intervals 


Time, 

sec 

h/k 

R, 

In. 

Reference 

[4] 

N = 10 
m = 4 
ft = 1 .0 

N = 19 
m = 8 
ft = 1.0 

5 

CO 

18 

0.77220 

0.773094 

0.772922 



10 

.01449 

.011353 

.011951 

10 

00 

18 

.84661 

.846719 

.846811 



10 

.11595 

.112112 

.113523 

30 

CO 

18 

.93546 

.935278 

.935521 



10 

.57722 

.575979 

.578198 

90 

CD 

18 

.99370 

.993686 

.993776 



10 

.95872 

.958596 

.959245 


TABLE IV. - COMPARISON OF EXPONENTIAL FINITE DIFFERENCE 
METHOD IN TWO-DIMENSIONAL CARTESIAN COORDINATES 
TO THE ALTERNATING DIRECTION IMPLICIT METHOD [4]. 

[For comparison to results In ref. [4] at x » y * o.] 


Time , 
sec 

Exact 

ADI 

[4] 

ft 3 1.0, N a 11, 
At » 0.01, m = 4 

Q - 1.0, N » 21, 
At = 0.0025, m . 9 

0.1 

0.09883 

0.09333 

0.09829 

0.09924 

.2 

.40354 

.40354 

.40256 

.40354 

.3 

.63179 

.63224 

.63080 

.63166 

.4 

.77486 

.77532 

.77403 

.77472 

.5 

.86252 

.86283 

.86187 

.86240 

.6 

.91607 

.91624 

.91559 

.91597 

.7 

.7 

.9487/ 

.94886 

.94841 
(Total of 
70 steps) 

—^—r = 0.2 
m ♦ 1 

.94869 
(Total of 
280 steps) 

— ^-7 » 0.1 
m + 1 
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TABLE V. - STEADY STATE HEAT TRANSFER IN TWO-DIMENSIONS 


Comparison of exponential finite difference technique to a Gauss-Seldel 
technique for the solution of Laplace's equation. 


X 

y 

9 by 9 

Grid 

5 by 5 Grid 

Gauss-Seldel 
88 Iteratlons(a) 

Exponential 
finite 
difference 
40 Iterations 

Gauss-Seldel 
22 Iteratlons(a) 

Exponential 
finite 
difference 
20 Iterations 

0.000 

1.000 

0.0 

0.0 

0.0 

0.0 

.125 

.875 

1.7413 

1.7414 



.250 

.750 

6.8946 

6.8949 

7.1428 

7.1430 

.375 

.625 

15.0330 

15.0335 



.500 

.500 

24.9999 

25.0004 

25.0000 

25.0003 

.625 

.375 

34.9667 

34.9672 



.750 

.250 

43.1052 

43.1055 

42.8571 

42.8573 

.875 

.125 

48.2587 

48.2588 



1 .000 

0.000 

100.000 

100.000 

100.000 

100.000 


a From reference [4]. 


TABLE VI. - COMPARISON OF EXPONENTIAL, PURE-EXPLICIT, AND 
IMPLICIT FINITE DIFFERENCE METHODS FOR ONE-DIMENSIONAL, 
UNSTEADY-STATE HEAT TRANSFER WITH TEMPERATURE VARYING 
THERMAL CONDUCTIVITY AT THE CENTER OF THE SLAB 


[K(T) = 1.0 * B(T) ; B > 0.01.] 


Time, 

sec 


Temperature, °C 


Exponential finite 
difference, N = 11 
m - 4, Q = 0.5, 
At = 0.005 sec 

Pure explicit, 

N = 11, Q = 0.25, 
At = 0.0025 sec 

Implicit, n . 1 .0 
At « 0.01 sec 

0.01 

98.15998 

100.00000 

94.35768 

.02 

88.87177 

89.21321 

85.90591 

.05 

61.30161 

60.09306 

61 .31385 

.1 

34.37147 

33.41929 

35.37178 
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TABLE VII. - COMPARISON OF THREE DIFFERENT, THREE-DIMENSIONAL UNSTEADY STATE 

HEAT TRANSFER SOLUTIONS 


[T(x,y, z,o) - 1.0; T(x,y,L,t) - T(x,L,z,t) * T(L,y,z,t) = 0; 8T/dx( o , y , z , t ) = 
dT/3y( x,o,z,t) = dT/az( x,y ,o, t) * 0; N * number of nodes In x, y, and z 
directions; ft = a At/(Ax)2 and Ax = Ay * Az.] 


Elapsed 

time, 

sec 

Position from 
center along 
diagonal 
x = y =* z 

Exact 

analysis 

result, 

°C 

Exponential finite 
difference results, 
°C 

N * 11, m * 4, 
ft = 0.75 

(a) 

Pure explicit 
finite difference 
results, 

°C 

ft a 0.15, N ■ 11 
(a) 

Method of Douglas 
Douglas finite 
difference results, 
°C 

ft a 0.15, Noll 

(a) 

0.09 

0.0 

0.893490 

0.892237 (0.14) 

0.889437 (0.45) 

0.886760 (0.75) 


.5 

.440712 

.440650 ( .014) 

.435058 (1.28) 

.439665 ( .24) 


.9 

.006491 

.006484 ( .11) 

.006319 (2.65) 

.006510 (-.29) 

.15 

0.0 

.645469 

.645209 ( .04) 

.640025 ( .84) 

.641484 ( .62) 


.5 

.253065 

.253286 (-.09) 

.250102 (1.17) 

.252641 ( .17) 


.9 

.003015 

.003022 (-.23) 

.002970 (1.49) 

.003023 (-.27) 


a Accuracy percent. 


TABLE VIII. - COMPARISON OF C.P.U. TIME ON TWO 
DIFFERENT MAINFRAMES FOR THREE DIFFERENT 
THREE-DIMENSIONAL FINITE DIFFERENCE METHODS 


[One-hundred time steps for each method.] 


Computer 

Exponential 8 

Method of 

Pure-explicit 


method , 

Douglas , 

method , 


sec 

sec 

sec 

CRAY-XMP 

0.2778 

0.955 

0.0627 

IBM-3033 

5.4 

12.6 

1 .8 


a Based on the number of sub-time Intervals 
equal to 100. 
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TABLE IX. - COMPARISON OF EXPONENTIAL 
FINITE DIFFERENCE METHOD TO EXACT 
RESULTS OF BOUNDARY LAYER EQUATION 
[TO] FOR THE VELOCITY PROFILE AT 
ONE DOWNSTREAM LOCATION. 

[Distance downstream x = 500 cm, 

V = 0.0072 cn)2/s.] 


Distance 
perpendicular 
to plate, 
y (cm) 

Exact 

result 

[10] 

Exponential 
method result, 
N = 21 m * 8 

1 .0 

0.17 

0.17428 

2.0 

.34 

.34643 

3.0 

.51 

.51020 

4.0 

.65 

.65658 

5.0 

.78 

.77684 

6.0 

.87 

.86636 

7.0 

.93 

.92638 

8.0 

.96 

.96265 


52 


i,j - SPACIAL GRID 
n - Tine GRID 



FIGURI 1.- COMPUIATIONAL GRID I OR IXPONtNTIAl F INI H DIFFt Rf NCF 
TECHNIQUt FOR 2-DlflENSIONAl CARTLSIAN COORDINAItS. 
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1st 2nd 

TIME SUB TIME SUB- 
INTERVAL INTERVAL 


n 

NODE i-1 • 


NODE i • 


NODE i+1 • 


TIME STEP 

n+l/3 n+2/3 n+1 

• T 

AX 









TOlAl DIWNSIONII SS 1 1 ME STEP 

FIGURE 2. COMPUTAI IONAI GRID FOR 2 IIMI SUB INII RVAI S <M - 2). CARI 
TESIAN COORDINATE. (SHOWN FOR 1 -DIME NSION) . 
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FIGURE 3. - FLOW DIAGRAM FOR 2-DIMENSIONAL EXPONENTIAL 
FINITE DIFFERENCE ALGORITHM. 

NTOT - TOTAL NUMBER OF TIME STEPS 
P, : - SUM OF DIMENSIONLESS DRIVE NUMBERS 
VTj'*:- FIELD VARIABLE DURING TIME STEP 
M, j - DIMENSIONLESS DRIVE NUMBER 
NSV - NUMBER OF TIME SUB- INTERVALS 
q A t 

n - DIMENSIONLESS TIME ( x) 2 
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FIGURE i*. - EFFECT OF NONDIMENSIONAL TIME STEP SIZE ON 1 NODE MODEL 
SOLUTION. (FROM REF. [2D. 


‘j6 




I IGURf S. T XHONl NI I Al FINITL DIFF FRF NCF NODAI LONF IGURAT ION FOR 
1 DIMFNSIONAI HI AT TRANSIT R WITH FINITE MEAT IRANSFf R COEFFICIENT. 



dT 

dr (R 1' t) = 0 

R 1 10.0 in. R ? = is.o in. 

I KiURr b. PROBUM CONDITIONS I OR COMPARISON Of IXPONtNTIAl 
flNINH DM 1 1 Rt NCI IICHNIOUl 10 CHARACII RISTIC PROBUM SOI UT ION. 


!>8 


I.C. : l(x,y,0) 0.0 



I IGURI /. PROBUM DtSCRIPT ION FOR 7 I) I ft NS I ONAl COMPARISON OF 
LXPONENIIAL I INI ft DIFFERENCE METHOD TO THE ALTERNATING DIRECTION 
IMPLICIT (ADI) METHOD. 
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T = 100 

FIGURE 8. PROBLEM SKLICH FOR COMPARISON OF EXPONTENTIAL FINITE 
DIFFERENCE SOLUTION TO THE GAUSS-StIDEL ITERATIVE METHOD FOR 
THE SOLUTION OF LAPLACES EQUATION IN 7-DIMENSIONS. 
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Kx.o) 100.0 

1 ( 0 , 1 ) 1 ( 1 . 1 ) 0.0 

K(f) K r ♦ 0K R T 


(A) ONL DIMKNSIONAl PROBLEM WITH VARYING THERMAL CON- 
DUCTIVITY. 





*R 

ILMPLRATURL 

(IT) I INIAR RHAIIONSHIP Itt IWI I N CONDUCTIVITY AND TtMPERATURt . 

I IGIIRI ‘1. !>KI I CHI S SHOWING PROITII M SIAM Ml NT FOR TEMPERATURl VARY- 
ING IIIIRMAI CONDUCT IV1 IY. 
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100 


80 


i 


6 



a 



c 


60 



20 


a 


A 

O EXPONENTIAL S 

□ PURE EXPLICIT 
A IMPLICIT 


& 

0 


J I I I i 

.2 .M .6 .8 1.0 

DIMENSIONLESS POSITION, x/L 


H IGURE 10. COMPAR I SON OF ME I HODS FOR TlMPERATURE-VARYING 
CONDUCTIVITY, SHOWING TtMPERATURE FIELD AT t = 0.02 SEC. 
K(T) K R (UPl), WHERE K R 1.0; p 0.01; T ( X ,0 ) = 100; 
1(0, t) KUt) = 0; a = 1. 
i XPONLNTIAL : M H; tt = 0.5; 20 TIME STEPS, 

PURI: IMPl IC1T: 12 = 0.25; 8 TIME STEPS. 


MOVING WALL VELOCITY = U, ILMPERATURt = T Q 



BOUNDARY CONDITIONS: T( y ,o) = T 0 x <0 

T< 0 ,x) = T-j X >0 

FIGURE 11. - SKETCH SHOWING CONDITIONS FOR DEVELOPING TEMPERATURE 
FIELD IN LAMINAR COUETTE FLOW. 
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FIGURE 1?. DFVtLOPING TEMPI RATURE FIFED IN IAMINAR COUETTE 
I I OW, SHOWING I HE EFFFCT OF FIUID VELOCITY "U"; (L * 1.0 m). 






3-DIMENSIONAL UNSTEADY STATE CONDUCTION 

fy 


/ ( 0 . 1 , 1 ) 


( 1 . 1 . 1 ) 



BOUNDARY CONDITIONS: t ^0 


dl dT dT 

dx <o,y.z.t) dy (x.o.z.t) - dz <x.y,o,t> - 0 


Td.y.Z.t) T(x,1,z,t) = T(x,y,1, t) = 0 
FIGURr 13. - BOUNDARY AND INITIAL CONDITIONS FOR THREE DIMENSIONAL 
UNSILADY STATE CONDUCTION HEAT TRANSFER. 
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FIGURE 1M. METHOD OF DOUGLAS SHOWING THE PROCEDURE USED TO SWEEP IN 
THE SUCCESSIVE COORDINATE DIRECTIONS. 


O EXPONENTIAL FINITE 

DIFFERENT RESUl T 
EXACT ANAIYSIS 



.? .4 .6 .8 1.0 


X/L 

I IGURt IV - COMPARISON OF SItADY STATE SOLUTIONS COMPARING 
I Hi IXACI Rl SUITS TO I ML EXPONENTIAL FINITE DIFFERENCE 
SOI III ION U(o, t ) 1.0; U(L,t) 0.0. 
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I KilJRt 16. iXPONINIAl FINI1I DIIFtRtNCC RtSULIS FOR VARY 
INO KlNtNMlC VISCOSITY. ALl VHOCmfcS ARE SHOWN FOR 

t 1.0 seconds WITH; N = 21. M * 8, — 7 = *1.0. 

(AXr 

U(o.t) * 1.0 

u(l. t) = o.o 


y 



OintR COND1 1 1ONS: Wx.o) 0; Vlx, 0 > 0 

Vlx.l ) 0 


I MilIKt 1/. IUHINDARY I AYl R IX VI l OPH N! Al ONI, A COOIFD II A I PI All . 


VELOCITY BOUNDARY LAYER THICKNESS, CM 
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riGURfc 18. - EXPONENTIAL FINIIE DIFFtRENCE RESULTS FOR BOUNDARY LAYtR tQUATIONS. WITH CONDITIONS U(o,y) = 1.0, 
IJ(x.o) 0, V(x,o) * 0. V(x,L ) 0, T(x, 0 ) = 0.0, T(0,y) = 1.0; V = 0.0072 cm 2 /sec; a = 0.01 cm 2 /sec. 


APPENDIX 


This appendix contains all the computer programs mentioned In this 
report. A computer program variable list Is also contained with a 
description of their use, and a program number to refer to the programs 
that they are contained In. 

Each of these programs was written to be run In an Interactive 
mode with the mainframe computer. The only cases run differently were 
for the three-dimensional unsteady state heat transfer cases that were 
run In batch mode on the Cray X-MP. 

The program structure Is as follows. A main program Is used to 
describe the necessary parameters and for asserting the proper boundary 
conditions. The main program then calls the subroutine where the 
actual finite difference methods are exercised and the results are then 
printed. 
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Program 

number 

1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 


Program 

name 

SOURCE. EFOCYL 

SOURCE. EFD20 

SOURCE. LAPLAC 
SOURCE. EFDVAR 

SOURCE. EXPVAR 

SOURCE. IMPVAR 

SOURCE. COUE 
SOURCE. EX3D 
SOURCE. EFD3D 

SOURCE. EXPL3D 

SOURCE. DOUGLA 

SOURCE. BURGER 
SOURCE. EXBURG 
SOURCE. NONBOU 


COMPUTER PROGRAM LIST 

Program function 


One-dimension, unsteady state, cylindrical 
coordinates. Infinite and finite heat 
transfer coefficient 

Two-dimensional Cartesian coordinates, 
unsteady state heat transfer 

Two-dimensional Laplace's equation 

One- dimensional unsteady state heat 

conduction, varying thermal conductivity, 
exponential finite difference method 

One-dimensional, unsteady state heat 
conduction, varying thermal conductivity, 
explicit finite difference method 

One-dimensional, unsteady state heat 
transfer, varying thermal conductivity. 
Implicit finite difference method 

One-dimensional, developing temperature field 
In laminar couette flow 

Exact analysis, three-dimensional heat 
transfer In a cube 


Three-dimensional unsteady state heat 

transfer In a cube using exponential finite 
difference method 


Three-dimensional unsteady state heat 
transfer In a cube using explicit finite 
difference method 


Three-dimensional unsteady state heat 
transfer In a cube using the method of 
Douglas 


Exponential solution of nonlinear viscous 
Burger's equation 


Pure explicit solution of nonlinear viscous 
Burger's equation 

Exponential Method of solution for boundary 
layer equations for flow over a flat plate 
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Program 

variable 

name 

N 

NS 

NTOT 

TS1 

T 

0L 

R 

I PR 
NB 
V 
HK 


VM 


M 

P 

VT 

B 

THE 

TIME 


UMAX 

ET1 


COMPUTER PROGRAM VARIABLE LIST 


Programs Variable Description 

used In 


1-14 

1-12,14 

1-14 

1-13 


Number of nodes 

Number of time sub Intervals 

Total number of time or spaclal steps 

Dimensionless time step Increment 


1-7,9,10,12,13 


1 


1 

1-14 

1 


1-14 


1 


Total elapsed time or spatial distance 
between steps 

Radial distance between adjacent nodes 
Radial length 

Number of steps between output of results 
Heat transfer boundary condition flag 
Dependent variable 

Convection heat transfer coefficient 
divided by thermal conductivity 


1 

1 - 10,12 

1 - 10,12 

1-14 

1 


Dependent variable value outside of 
solid In the surrounding medium 

Dimensionless drive number 

Sum of the drive numbers 

Dependent variable value during 
sub-time Interval 

Blot number 


1-3,12-14 Variable used for the output of results 

1-6,8-13 Total elapsed time of the solution at 

the current output 


1-14 Output counter 

3 Accuracy desired In solution of 

Laplace's equation 
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Program 

variable 

name 

Program 
numbers 
used In 

Description 

DELV 

3 

Difference In dependent variable value from 
one time step to the next 

ETA 

3 

Sum of the absolute value of the differences 
found In DELV 

KR 

4-6 

Reference thermal conductivity 

BETA 

4-6 

Slope of thermal conductivity variation with 
temperature 

KO , K 

4 

Thermal conductivity at the total time step 
Interval or sub-time Interval respectively 

THE 

4 

Klrchoff transformation variable (used In 
exponential finite difference program with 
varying thermal conductivity) 

KAPPA 

4 

All values known from the last time step 
Increment and used to solve the quadratic 
equation that results In the exponential 
finite difference solution with temperature 
varying thermal conductivity 

DERSQR 

5 

Absolute value of velocity difference found In 
evaluation of velocity gradient 

OMEGA 

6 

Same as nondlmenslonal time step 

A,B,C,D 

6,11 

Coefficient used In trldlagonal matrix 
algorithm. 

KAP.GAM 

6 

Variables used to determine A,B,C 

BETA, GAMMA 

6 

Variables used In Thomas, trldlagonal 
algorithm 

TS 

6 

Same as total elapsed time 

T 

6 

Dependent variable 

V 

6 

Solution vector tri-diagonal algorithm 

SP 

7 

Maximum width divided by maximum velocity 

FL 

7 

Parameter based on position In flow 

DIST 

7 

Serves same function as time for unsteady 
state problem 
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Program 

variable 

name 

Program 
numbers 
used In 

Description 

PI , PI2, PI3 

8 

ir, ir^, 

NOOES 

8 

Same an N 

TO 

8 

Initial temperature 

TI 

8 

Surface temperature, t > 0 

ALFA 

8 

Thermal dlffuslvlty 

TNEW 

8 

Exact temperature at a x, y, and z 

location after elapsed time t has occurred 

T 

11 

Dependent variable 

DELX.DELY.OELZ 

11 

Part of the central difference operator 

u,v 

11 

Variables used to sweep preliminary solution 
In the x then y directions respectively 

M 

11 

Used as an array to contain the known 
quantities used In the trldlagonal algorithm 

RNU 

12-14 

Kinematic viscosity 

DX 

14 

Step length In flat plate direction 

RAL 

14 

Thermal dlffuslvlty 

YMAX 

14 

Maximum distance perpendicular to flat plate 

DY 

14 

Step length perpendicular to the flat plate 

u,v 

14 

Dependent variables (velocity) In x and y 
directions respectively 

T 

14 

Temperature field variable 

MU.MT 

14 

Drive numbers for velocity In x-dlrectlon 
and temperature field 

PU.PT 

14 

Sum of drive numbers for MU, MT 

UT,TT, VT 

14 

x-dlrectlon velocity, temperature, and 
y-dlrectlon velocity on subinterval 

UT1 

14 

Temporary U-dlrectlon velocity field 

TS.TS1 

14 

Dimensionless time step for temperature and 
x-dlrectlon velocity respectively 


non non ooonoooooooooo 
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WRITTEN BY R.F. HANDSCHUH 

SOURCE. EFDCYL 

***** PROGRAM #1 ***** 


THIS PROGRAM IS TO BE USED AS THE STARTING POINT FOR INVESTIGATING 
THE EXPONENTIAL FINITE DIFFERENCE ALGORYTHM. THIS METHOD WAS INTRODUCED 
BY M.C. BHATTACHARYA. THIS SOURCE CODE IS FOR CYLINDRICAL 
COORDINATES, UNSTEADY-STATE HEAT CONDUCTION, 1 DIMENSION. 


IMPLICIT REAL*8(A-H,0-Z) 

REAL*8 V(IOO) , RC 1 0 0 ) 

INPUT THE PROGRAM DATA 

WRITEC 6,15) 

15 FORMAT( IX, 'NUMBER OF NODES=N 13'/) 

READ( 9 , 1 0 )N 

10 FORMATCI3 ) 

WRITEC6 , 12) 

12 FORMATC IX, 'NUMBER OF TIME SUB INTERVALS= NS 13') 

P.EAD( 6 , 13 )NS 

13 ’ FORMATCI3) 

URITEC6 , 16 ) 

16 FORMATC IX , ’ TOTAL NUMBER OF TIME STEPS= NTOT 13') 

READ( 9,21 )NTOT 

21 FORMATC 13 ) 

WRITEC 6,24) 

24 FORMATC IX, ’INPUT TIME*THERMAL DIFFUSIVITY / RAD INT SQUARED F5.3') 
READC 9 , 25 )TSI 

25 FORMAT CF5 .3) 

WRITEC 6, 22) 

22 FORMATC IX, 'TOTAL TIME OF ONE TIME STEP= T F7 . 5 ' ) 

READC 9, 23 )T 

23 FORMAT CF7. 4) 

WRITEC 6, 26) 

26 FORMATC IX, ’INPUT RADIAL INTERVAL I.ENGTH=DL F5.3') 

READ C 9 , 27 ) DL 

27 FORMATC F5. 3) 

RC 1 ) = 1 . 

DO 28 1=2, N 
IM1=I-1 

28 RCI)=RCIM1)-DL 
WRITEC 6, 32) 

32 FORMATC IX, 'INPUT NUMBER OF TIME STEPS BEFORE PRINTING 13') 

P.EADC 9,33 )IPR 

33 FORMATCI3) 

DETERMINE THE TYPE OF BOUNDRY CONDITION, THEN SET VALUES 
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WRITEC 6,14) 

14 FORMATC IX, 'INPUT HEAT TRANSFER B.C. 0 - INFINITE 1 - FINITE'/) 
READ( 9 , 17 )NB 

17 F0RMATCI1) 

IFCNB.EQ.l) GO TO 100 
V(1)=0. 

DO 30 1=2, N 

30 V(I)=1.0 

CALL EXP FIN DIF FOR INFINITE HEAT TRANSFER COEFFICENT 

CALL EFDIHC ( N , NS , NTOT , TSI , V , T , R , DL , IPR ) 

GO TO 101 

100 CONTINUE 
WRITEC 6, 31) 

31 FORMAT ( IX, 'INPUT HEAT COEF / TERM COND F5.3') 

READC 9 , 34 )HK 

34 FORMAT ( F5 . 3 ) 

DO 40 1=1, N 
40 V(I)=1. 

VM=0 . 

CALL EXP FIN DIF FOR FINITE HEAT TRANSFER COEFFICENT 

CALL EFDFHC ( HK , N , VM , DL , NS , V , NTOT , TSI , T , R , IPR ) 

101 CONTINUE 
STOP 
END 

SUBROUTINE EFDIHC 

SUBROUTINE EFDIHC ( N , NS , NTOT , TSI , V , T , R , DL , IPR ) 

FOR INFINITE HEAT TRANSFER COEFFICENT 
IMPLICIT REAL*8 ( A-H , O-Z) 

REAL*8 VT(100),V(100),M(100),P(100),R(100), THE (100) 
TS=TSI/DFLOAT(NS+l ) 

WRITE (6,21)(R(I),I=1,N) 

21 FORMAT ( 1X,11(F6.3,2X)) 

WRITEC 6,22 )DL , TSI 

22 FORMAT ( 1X,2(F6.3,2X)) 

N1=N-1 

NS1=NS+1 

BEGIN MAIN TIME STEP LOOP 
DO 20 L=1 , NTOT 

ZERO DRIVE NUMBERS AND SET TEMPORARY VARIABLES EQUAL TO THE 
LAST TOTAL TIME STEP VALUES 

DO 15 1=1, N 

15 P(I)=0. 


non ooo non non ooo ooo 


75 


DO 10 1=1, N 
10 VT(I)=V(I) 

SUB TIME INTERVAL 

DO 30 K=1 , NS1 

CALCULATE THE DRIVE NUMBERS 

DO 40 1=2, N1 

IM1=I-1 

IP1=I+1 

M( I ) = ( 2 . *VT CD-VTCIMl )-VTCIPl ) ) /VT C I ) 

40 M(I)=M(I)+DL*( VTCIPl )-VT(IMl ) )/( VT(I)*2 .*RCI) ) 

M(N)=2.*CVT(N)-VT(N1) )/VT(N) 

CALCULATE THE DEPENDENT VARIABLE ON THE SUB-INTERVAL LEVEL 
DO 50 11=2, N 

50 VTCI1)=VTCI1)*DEXPC-TS*MCI1)) 

SUM THE DRIVE NUMBERS 

DO 60 1=2, N 
60 P(I)=P(I)+MCI) 

30 CONTINUE 

CALCULATE THE DEPENDENT VARIABLE ON THE NEXT COMPLETE TIME STEP 
DO 7 0 1=1, N 

V ( I ) =V( I ) *DEXP ( -TS*PCI ) ) 

70 CONTINUE 
ITMAX=ITMAX+1 

PRINT THE RESULTS 

IF ( ITMAX . LT . IPR ) GO TO 20 
DO 71 1=1, N 

71 THE(I)=V(I) 

ITMAX=0 
WRITEC 6,5) 

5 FORMAT (/ ) 

WRITE(6,31)L 

31 FORMAT( IX, 'TIME STEP NUMBER=',I3) 

TIME=T*DFLOAT CL) 

WRITEC 6 , 32 ) TIME 

32 FORMATC IX, ’ELAPSED TIME= ' , FI 0.4, ' SECONDS') 

IFCN.GT. 11 )N21=N/2 

IFCN.GT. 11 )GO TO 81 
WRITE C6,82)( THE (I),I=1,N) 

82 FORMATC 11(2X,F8.6)) 

GO TO 84 
81 CONTINUE 

WRITEC 6,82) CTHECI) ,1=1 ,N21 ) 


lb 


NNS=N2 1 + 1 

WRITE ( 6 , 82 ) ( THE ( I ) , I=NNS , N ) 

84 CONTINUE 
20 CONTINUE 
RETURN 
END 
C 

C SUBROUTINE EFDFHC 

C 

SUBROUTINE EFDFHC ( HK , N , VM , DL , NS , V , NTOT , TSI , T , R , IPR ) 

IMPLICIT REAL*8 (A-H,0-Z) 

REAL*8 VT(100),V(100),M(100),P(100),R(100), THE (100) 

B=HK*DL 

TS=TSI/DFLOAT( NS+ 1 ) 

NS1 =NS+ 1 
N 1 = N- 1 
N2=N“2 
C 

C BEGIN THE TOTAL TIME STEP LOOP 

C 

DO 20 L=l,NTOT 
C 

C ZERO THE DRIVE NUMBERS AND SET THE TEMPOARY DEPENDET VARIABLES 

C EQUAL TO THE LAST COMPLETE STEP VALUES 

C 

DO 15 1=1, N 
15 P(I)=0. 

DO 10 1=1, N 
10 VT ( I ) =V ( I ) 

C 

C SUB TIME INTERVAL 

C 

DO 30 K=1 , NS1 
C 

C CALCULATE THE DRIVE NUMBERS 

C 

M(l)=~(2. *VT (2)-(2.+2.*B)*VT(l)+2.*B*VM)/VT(l) 

M( 1 )=M( 1 )-DL^B^(VT( 1 )-VM)/(R( 1 )*VT( 1 ) ) 

DO 40 1=2, N1 
IM1 =I~ 1 
IP1=I+1 

M(I)=(2.*VT(I)-VT(IM1)-VT(IP1) )/VT(I) 

40 M(I)=M(I) +DL*( VT(IP1 ) -VT( IM1 ))/(VT(I)^2.^R(I) ) 

M(N)=2.*(VT(N)-VT(N1))/VT(N) 

C 

C CALCULATE THE DEPENDENT VARIABLE ON THE SUB-INTERVAL LEVEL 
C 

DO 50 11=1, N 

50 VT(I1 )=VT(I1)*DEXP(-TS*M(I1) ) 

C 

C SUM THE DRIVE NUMBERS 
C 

DO 60 1=1, N 
60 P(I)=P(I)+M(I) 


non ooo 
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30 CONTINUE 

CALCULATE THE DEPENDENT VARIABLE ON THE NEXT COMPLETE TIME STEP 
DO 70 1=1, N 

V ( I ) =VC I )*DEXP( -TS*PCI ) ) 

70 CONTINUE 
ITMAX=ITMAX+1 
IF(ITMAX.LT. IPR)GO TO 20 

PRINT THE RESULTS 

DO 71 1=1, N 

71 THE(I)=V(I) 

ITMAX=0 
WRITEC 6,5) 

5 FORMATC/ ) 

WRITEC 6 , 3 1 ) L 

31 FORMATC IX, ’TIME STEP NUMBERS, 13) 

TIME=T*DFLOATC L ) 

WRITEC 6 , 32 ) TIME 

32 FORMATC IX , ' ELAPSED TIME= ' , FI 0 . 4 , ’ SECONDS’) 

IFCN.GT.il )N21=N/2 

IFCN.GT. 11 )GO TO 81 
WRITE C 6 , 8 2 ) C THE C I ) , 1= 1 , N ) 

82 FORMATC 11C2X,F8.6)) 

GO TO 84 
81 CONTINUE 

WRITEC 6 , 82 ) C THEC I ) , 1=1 , N21 ) 

NNS=N2 1+1 

WRITEC 6 ,82) C THEC I) ,I=NNS,N) 

84 CONTINUE 
20 CONTINUE 
RETURN 
END 
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SOURCE. EFD2D 

WRITTEN BY R.F. HANDSCHUH 
***** PROGRAM #2 ***** 


THIS PROGRAM IS FOR 2-DIMENSIONAL CARTESIAN COORDINATES 
UNSTEADY STATE HEAT TRANSFER. THE METHOD OF SOLUTION IS THE 
EXPONENTIAL FINITE DIFFERENCE ALGOR YTHM. THIS PARTICULAR PROGRAM IS 
FOR INFINITE HEAT TRANSFER COEFFICENT AT THE EXPOSED SURFACES 
AT X=Y= 1 . 0 FOR X=Y=0 THE SURFACE IS CONSIDERED TO BE 
PERFECTLY INSULATED. 


IMPLICIT REAL*8(A-H,0-Z) 
RF.AL*8 VC 25, 25) 

INPUT PROGRAM DATA 


5 

0 


13 

16 

21 

24 

25 
22 
23 

26 
27 
250 


251 


WRITE C 6 ,15) 

FORMATC IX, 'NUMBER OF NODES=N 13'/) 

READC9, 10 )N 
FORMAT (13) 

WRITEC 6,12) 

FORMATC IX, 'NUMBER OF TIME SUB INTERVALS 2 NS 13') 

READC 9 , 13 )NS 
FORMATCI3 ) 

WRITEC6 , 16 ) 

FORMATC IX, 'TOTAL NUMBER OF TIME STEPS 2 NTOT 13') 

READC 9,21 )NTOT 
FORMATCI3) 

WRITEC 6, 24) 

FORMATC IX, 'INPUT TIME*THERMAL DIFFUSIVITY / LENGTH SQUARED 
READC 9 , 25 )TSI 
FORMATC F5 . 3 ) 

WRITEC 6 ,22) 

FORMATC IX, 'TOTAL TIME OF ONE TIME STEP 2 T F6.4’) 

READ ( 9 , 23 )T 
FORMATC F6. 4) 

WRITEC 6 ,26 ) 

FORMATC IX, 'NUMBER OF TIME STEPS BEFORE PRINTING 2 13') 

READC 9,27 )IPR 
FORMAT (13) 

WRITEC 6 ,250 )N, NS, NTOT 

FORMATC IX,'# OF NODES 2 ', 13 , 2X, ' # OF SUB-TIME-INT 2 ' , 13 , 2X , 
*'# OF TIME STEPS 2 ',13) 

WRITEC 6,251 )TSI , T 

FORMATC IX, ' (TIME*THER DIFF) /LENGTH SQUARED 2 ' F5 . 3 , 2X , 

* ' TIME STEP LENGTH 2 ' ,F6 .4/) 

N1=N-1 


F5 . 3 ' ) 


INITIALIZE THE BOUNDRY CONDITIONS 
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DO 30 1=1, N1 
DO 30 J=1,N1 
30 V(I,J)=1.0 

I=N 

DO 50 J=1,N 

50 V(I,J)=0.0 
J=N 

DO 51 1=1, N 

51 V(I,J)=0.0 

CALL EXP FIN DIF FOR INFINITE HEAT TRANSFER COEFFICENT 

CALL EFDIHC ( N , NS , NTOT , TSI , V , T , IPR ) 

STOP 

END 

SUBROUTINE EFDIHC 

SUBROUTINE EFDIHC ( N , NS , NTOT , TSI , V, T, IPR ) 

FOR INFINITE HEAT TRANSFER COEFFICENT 
IMPLICIT REAL*8 (A-H,0~Z) 

REAL*8 VT(25,25) , VC 25,25 ) ,M( 25,25 ) ,P( 25,25 ) ,THE( 25, 25) 

PRINT HEADING 
WRITEC 6,222) 

222 FORMAT (IX , ************ SOURCE. EFD2D a******#****'//) 

TS=TSI/DFLOAT(NS+l) 

N1=N-1 

NS1=NS+1 

BEGIN MAIN TIME STEP LOOP 
DO 20 L=1 , NTOT 

ZERO THE DRIVE NUMBERS AND SET TEMPORARY DEPENDENT VARIABLE 
EQUAL TO THE LAST FULL TIME STEP VALUE 

DO 15 J=1 , N 
DO 15 1=1, N 
15 P(I,J)=0. 

DO 10 J=1 , N 
DO 10 1=1, N 
10 VT(I,J)=V(I,J) 

SUB TIME INTERVAL 

DO 30 K=1 ,NS1 

CALCULATE THE DRIVE NUMBERS 
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DO 41 J=2,H1 

JM1=J-1 

JP1=J+1 

DO 40 1=2 , N1 

IM1=I-1 

IP1=I+1 

MCI, J)=(VTCIP1, JJ+VTCIM1, JJ+VTCI, JP1J+VTCI, JM1 )-4 .*VT(T, J) ) 

40 MCI, J)=M(I, J)/VT(I, J) 

41 CONTINUE 

INSULATED BOUNDRY ALONG X-AXIS 
J = 1 

JP1=J+1 
DO 42 1=2, N1 
IP1*I+1 
IM1 =1- 1 

42 MCI, J)=(VT(IP1 , J)+VT(IM1, J)+2.*VT(I, JP1 ) -4 . *VT( I , J) )/VT(I, J) 
INSULATED BOUNDRY ALONG Y-AXIS 

1=1 

IP1=I+1 
DO 43 J=2,N1 
JP1=J+1 
JM1 = J- 1 

43 MCI, J)=C2.*VTCIP1, JJ+VTCI, JP1J+VTCI, JM1 ) -4 . *VTCI, J))/VT(I, J) 
CORNER AT ORIGIN 

M ( 1 , 1 ) = C 2 . *VT C 1 , 2 ) +2 . *VT C 2 , 1 ) -4 . *VT C 1 , 1 ) ) /VT C 1 , 1 ) 

CALCULATE THE DEPENDENT VARIABLE ON THE SUB-INTERVAL LEVEL 

DO 50 11=1, N1 
DO 50 J1 =1 , N1 

50 VTCI1, J1)=VT(I1 , J1)*DEXP(TS*MCI1, Jl) ) 

SUM THE DRIVE NUMBERS 

DO 60 1=1, N1 
DO 60 J= 1 , N1 

60 PCI, J)=PCI, JJ+MCI, J) 

30 CONTINUE 

CALCULATE THE DEPENDENT VARIABLE AT THE NEXT COMPLETE TIME STEP 

DO 70 J=1 , N 
DO 70 1=1, N 

V(I, J)=V(I, J)*DEXP(TS*P(I, J) ) 

70 CONTINUE 

ITMAX=ITMAX+1 

IF ( UMAX . LT . IPR )GO TO 20 


on 


HI 


PRINT THE RESULTS 

WRITE (6,5) 

5 FORMATC/) 

WRITEC 6 , 31 )L 

31 FORMAT (IX, 1 TIME STEP NUMBER=M3/) 
TIME=T*DFLOATC L ) 

WRITE ( 6,32 )TIME 

32 FORMATC 5X, ’ELAPSED TIME= ' , FI 0 . 4 , 1 SECONDS ’ / ) 
DO 71 1=1, N 

DO 71 J=1,N 

71 THECI, J)=l . O-VCI, J) 

IFCN.GT. 11 )GO TO 58 
DO 59 J=1,N 

WRITEC 6, 82 ) (THECI, J), 1=1, N) 

82 FORMAT (11(2X,F8.6)) 

59 CONTINUE 
GO TO 54 
58 CONTINUE 

DO 57 J=1 , N 

WRITEC 6 , 5 6 ) C THE Cl, J), 1=1, 11) 

56 FORMAT C11(2X,F8.6)) 

WRITEC 6 , 56 ) ( THECI, J) , 1=12 , N ) 

57 CONTINUE 
54 CONTINUE 

ITMAX= 0 
20 CONTINUE 
RETURN 
END 
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SOURCE. LAPLAC 

WRITTEN BY R.F. HANDSCHUH 

***** PROGRAM #3 ***** 


THIS PROGRAM IS TO BE USED TO SOLVE THE LAPLACE’S EQUATION 
USING THE EXPONENTIAL FINITE DIFFERENCE METHOD. 


IMPLICIT REAL*8(A-H,0-Z) 

REAL*8 V( 25 > 25 ) 

INPUT PROGRAM DATA 

WMTEC6 , 15) 

15 FORMATC IX, 'NUMBER OF NODES=N 13'/) 

READ( 9, 10 )N 

10 FORMAT (13) 

WRITEC6 , 12) 

12 FORMATC IX, 'NUMBER OF TIME SUB INTERVALS 2 NS 13') 
READ( 9 , 1 3 )NS 

13 FORMAT (13) 

WRITEC 6,16) 

16 FORMATC IX, 'TOTAL NUMBER OF TIME STEPS 2 NTOT 13') 

• READ (9,21) NTOT 

21 FORMATCI3 ) 

WRITEC 6, 24) 

24 FORMATC IX, 'INPUT TIME / LENGTH SQUARED F5.3') 

READ (9,25) TSI 

25 FORMATCF5 . 3) 

WRITEC 6,22) 

22 FORMATC IX, 'TOTAL TIME OF ONE TIME STEP= T F7.6') 

READC 9 , 23 )T 

23 FORMAT (F7 . 6 ) 

WRITEC 6, 26) 

26 FORMATC IX, 'NUMBER OF TIME STEPS BEFORE PRINTING 2 13’) 
READ (9,27 )IPR 

27 FORMATCI3 ) 

WRITEC 6, 31) 

31 FORMATC IX, 'INPUT ACCURACY DESIRED 2 F7 . 6 ’ ) 

READC 9, 32)ET1 

32 FORMATC F7 . 6 ) 

N1=N-1 

INITIALIZE THE BOUNDRY CONDITIONS 

DO 30 1=2, N1 
DO 30 J=2 , N1 
30 V(I,J)=100. 

C 


I=N 

DO 50 J=1 , N 
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50 V(I,J)=0.0 
J=N 

DO 51 1=1, N 

51 V(I,J)=0.0 


J=1 

DO 53 1=2, N1 
53 V(I,J)=0.0 


1=1 

DO 52 J=1,N 
52 V(I,J)=100.0 


CALL EXP FIN DIF FOR LAPLACE EQUATION 

CALL EFDIHC ( N , NS , NTOT , TSI , V , T , IPR , ET 1 ) 

STOP 

END 

SUBROUTINE EFDIHC 

SUBROUTINE EFDIHC ( N , NS , NTOT , TSI , V , T , IPR , ET1 ) 


IMPLICIT REAL*8(A-H,0-Z) 

REAL*8 VT(25,25) ,V( 25,25 ) ,M( 25,25 ) , PC 25,25 ) ,THE( 25, 25) 
TS=TSI/DFLOAT ( NS+ 1 ) 

N1=N-1 
NS1 =NS+ 1 

BEGIN MAIN STEP INCREMENT 

DO 20 L=1 , NTOT 
ETA= 0 . 0 

ZERO THE SUM OF THE DRIVE NUMBERS 

DO 15 J=1,N 
DO 15 1=1, N 
15 P(I,J)=0. 

SAVE THE DEPENDENT VALUES FROM THE LAST TOTAL TIME STEP 

DO 10 J= 1 , N 
DO 10 1=1, N 
DELV=V ( I,J)-VT(I,J) 

ETA=ETA+DABS( DELV ) 

10 VT ( I , J ) =V ( I , J ) 

IFCL.LE.88) GO TO 107 
IF ( ETA . LE . ET1 )GO TO 100 
107 CONTINUE 
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C 

C SUB TIME INTERVAL 

C 

DO 30 K=1,NS1 
C 

C CALCULATE THE DRIVE NUMBERS 

C 

DO 41 J=2,N1 
JM1= J- 1 
JP1=J+1 
DO 40 1=2, N1 
IM1=I-1 
IP1=I+1 

MCI, J)=(VTCIP1, J)+VTCIM1, JJ+VTCI, JP1)+VT(I, JM1 ) -4 . *VTC I , J) ) 

40 M(I, J)=M(I, J)/VT(I, J) 

41 CONTINUE 
C 

C CALCULATE THE DEPENDENT VARIABLE ON THE SUB-TIME INTERVAL 
C 

DO 50 11=2, N1 
DO 50 J1=2,N1 

50 VT(I1 , J1)=VT(I1, J1 )*DEXPCTS*MCI1, Jl) ) 

C 

C SUM THE DRIVE NUMBERS 

C 

DO 60 1=1, N1 
DO 60 J=1,N1 

60 PCI, J)=PCI, JJ+MCI, J) 

30 CONTINUE 
C 

C CALCULATE THE DEPENDENT VARIABLE ON THE NEXT TOTAL STEP 
C 

DO 70 J=1,N 
DO 70 1=1, N 

V(I, J)=V(I, J)*DEXP(TS*P(I, J) ) 

70 CONTINUE 
ITMAX=ITMAX+ 1 

C 

C PRINT THE RESULTS 

C 

IF(ITMAX.LT. IPR)GO TO 20 
WRITEC 6,5) 

5 FORMATC / ) 

WRITEC 6 , 3 1 ) L 

31 FORMATC IX, 'TIME STEP NUMBER= 1 ,13/ ) 

TIME = T*DFLOATC L ) 

DO 71 1=1, M 
DO 71 J=1,N 

71 THECI, J)=VCI, J) 

IFCN.GT. 11 )GO TO 58 
DO 59 J=1,N 

WRITEC 6,82) (THECI, J) ,1=1 ,N) 

82 FORMAT (11C2X,F8.4)) 

59 CONTINUE 


8b 


GO TO 54 
58 CONTINUE 

DO 57 J=1,N 

WRITE ( 6,56) (THECI, J) ,1=1,11) 

56 FORMAT (11(2X,F8.4)) 

WRITEC6 / 56 ) CTHECI/ J ) ,1=12/ N) 

57 CONTINUE 
54 CONTINUE 

ITMAX= 0 
20 CONTINUE 
GO TO 101 

100 CONTINUE 
WRITEC6 , 103)L 

103 FORMAT ( 2X 1 ***** CONVERGED RESULT ***** * , 5X , 1 ITERATION= ' , 14// ) 
DO 102 J=1,N 

WRITE ( 6,56)(V(I,J) , 1=1 ,N) 

102 CONTINUE 

101 CONTINUE 
RETURN 
END 
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WRITTEN BY R.F. HANDSCHUH 
***** PROGRAM #4 ***** 


SOURCE. EFDVAR 


THIS PROGRAM IS FOR THE SOLUTION OF THE DIFFUSION EQUATION 
WITH VARYING THERMAL CONDUCTIVITY. THE METHOD USED IS THE 
EXPONENTIAL FINITE DIFFERENCE METHOD. 


IMPLICIT REAL*8( A-H,0-Z) 
REAL*8 V ( 1 0 0 ) , KR 

INPUT PROGRAM DATA 


15 
10 
12 

13 

16 
21 

24 

25 
22 
23 

14 
17 

26 
27 


WRITEC6 , 15) 

FORMAT ( IX, 'NUMBER OF NODES=N 13'/) 

READ( 9 , 10 )N 
FORMAT (13) 

WRITE(6, 12) 

FORMAT ( IX, 'NUMBER OF TIME SUB INTERVALS 3 NS 13’) 

READ( 6 , 1 3 )NS 
FORMAT (13) 

WRITE( 6,16) 

FORMAT ( IX, 'TOTAL NUMBER OF TIME STEPS 3 NTOT 13') 

READ( 9.2DNTOT 
FORMAT( 13 ) 

WRITE( 6,24) 

FORMAT ( IX, 'INPUT TIME*THERMAL DIFFUSIVITY / LENGTH SQUARED 
READ( 9,25 )TSI 
FORMAT (F5. 3) 

WRITE( 6,22) 

FORMAT (IX, 'TOTAL TIME OF ONE TIME STEP 3 T F7.5') 

READ( 6 , 23 )T 
FORMAT(F7 .5) 

WRITE ( 6 , 14 ) 

FORMAT( IX, 'INPUT REFERENCE THERMAL CONDUCTIVITY=F5 . 4 ' ) 
READ( 9 , 17 )KR 
FORMAT(F5 . 4 ) 

WRITE( 6,26) 

FORMAT( IX, 'INPUT THERMAL CONDUCTIVITY SLOPE VALUE 3 F5 . 4 ' ) 
READ( 9 , 27 )BETA 
FORMAT( F5 . 4 ) 


F5 . 3 ' ) 


INITIALIZE THE BOUNDRY CONDITIONS 


V( 1 ) = 0 . 

V(N) = 0 . 
N1=N-1 

DO 30 1=2, N1 
30 V(I)=100. 
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CALL EXP FIN DIF FOR INFINITE HEAT TRANSFER COEFFICENT 

CALL EFDIHC C N , NS , NTOT , TSI , V , T , KR , BETA ) 

STOP 

END 

SUBROUTINE EFDIHC 

SUBROUTINE EFDIHC ( N , NS , NTOT , TSI , V , T , KR , BETA ) 

FOR INFINITE HEAT TRANSFER COEFFICENT 

IMPLICIT REAL*8(A-H,0-Z) 

REAL*8 VTC100) , V( 100 ) ,M( 100 ) ,P( 100 ) 

REAL*8 KR , K ( 1 0 0 ) , KO ( 1 0 0 ) , THE (100), KAPPA (100) 
TS=TSI/DFLOAT(NS+l ) 

N1=N-1 

NS1=NS+1 

PRINT HEADING 

WRITE( 6,5) 

WRITE ( 6,100) 

100 FORMATdX, '***** SOURCE . EFDVAR *****'/) 

BEGIN MAIN TIME STEP LOOP 
DO 20 L=1 , NTOT 

ZERO THE SUM OF THE DRIVE NUMBERS 

DO 15 1=1, N 
15 P(I)=0. 

SET VARIABLES EQUAL TO THE LAST STEPS VALUES 

DO 10 1=1, N 

10 VT(I)=V(I) 

DO 11 1=1, N 

11 K0 (I ) =KR+BETA*KR*V (I) 

SUB TIME INTERVAL 

DO 30 KK=1 , NS1 

DO 35 1=1, N 

K(I) =KR+BETA*KR*VT( I) 

35 THE(I)=VT(I)+BETA^VT(I)^VT(I)/2. 0 

CALCULATE THE DRIVE NUMBERS 

DO AO 1=2, N1 

IM1=I-1 

IP1=I+1 
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M(I)=VT(IP1)+VT(IM1 )-2.*VT(I) 

M( I ) =M( I ) +BETA*( VTC IP1 )**2 . +VTC IM1 ) **2 . -2 . *VT( I )**2 . )/2 . 

40 M( I ) =M( I )/THE( I ) 

DO 50 11=2, HI 

50 KAPPA (II )=THE(I1 )*DEXP(TS*( 1 . +BETA*VT(I1 ) )*M(I1 ) ) 

C 

C CALCULATE THE DEPENDENT VARIABLE ON THE SUB-TIME INTERVAL 
C 

DO 55 11=2, N1 

55 VTCIl )=(-l . +SQRTC 1 .+2 .*KAPPA(I1 )*BETA) )/BETA 
C 

C SUM THE DRIVE NUMBERS 

C 

DO 60 1=2, N1 
60 P(I)=P(I)+M(I) 

30 CONTINUE 
WRITE( 6,5) 

5 FORMAT ( / ) 

WRITEC 6 , 31 )L 

31 FORMATC1X, ’TIME STEP NUMBER= ' , 13/ ) 

TIME=T*DFLOAT ( L ) 

WRITEC 6 , 32 )TIME 

32 FORMATC5X, 'ELAPSED TIME= FI 0 . 4 SECONDS V ) 

C 

C CALCULATE THE NEXT TOTAL STEP DEPENDENT VARIABLES AND PRINT RESULTS 
C 

DO 70 1=1, N 

THE(I)=V(I) +BETA#V (I)*V(I)/2. 0 

KAPPA ( I ) =THE ( I ) *DEXP ( TS* ( 1 . +BETA*V (I) )*P(I) ) 

70 V(I) = (-1 .+SQRTC1 .+2.*KAPPA(I)*BETA»/BETA 

WRITEC 6,81)(V(I),I=1,11) 

81 FORMAT (1X,11(F8.5,2X)) 

20 CONTINUE 
RETURN 
END 
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WRITTEN BY R.F. HANDSCHUH 


SOURCE. EXPVAR 

***** PROGRAM #5 ***** 


THIS PROGRAM IS FOR THE SOLUTION OF THE DIFFUSION EQUATION 
WITH VARYING THERMAL CONDUCTIVITY. THE METHOD USED IS THE 
EXPLICIT FINITE DIFFERENCE METHOD. 


IMPLICIT REAL*8(A-H,0-Z) 
REAL*8 V(100),KR 

INPUT PROGRAM DATA 


15 
10 

16 
21 

24 

25 
22 


14 

17 

26 

27 


WRITE(6 , 15) 

FORMAT ( IX, ’NUMBER OF NODES=N 13’/) 

READ( 9 , 1 0 )N 
FORMATCI3) 

WRITE( 6,16) 

FORMAT( IX, ’TOTAL NUMBER OF TIME STEPS= NTOT 13*1 

READ( 9,21 )NTOT 

FORMATCI3) 

WRITEC 6,24) 

FORMATC IX, ’INPUT TIME*THERMAL DIFFUSIVITY / LENGTH SQUARED 
READ( 9 , 25 )TSI 
FORMAT (F5. 3) 

WRITEC 6,22) 

FORMATC IX, ’TOTAL TIME OF ONE TIME STEP= T F7.5’) 

READC 6 , 23 )T 
FORMATC F7 .5) 

WRITEC6, 14) 

FORMATC IX, ’INPUT REFERENCE THERMAL CONDUCTIVITY=F5 . 4 * ) 
READC 9 , 1 7 )KR 
FORMATCF5 . 4 ) 

WRITEC 6 ,26 ) 

FORMATC IX, ’INPUT THERMAL CONDUCTIVITY SLOPE VALUE= F5 . 4 ’ ) 
READ (9, 27) BETA 
FORMAT CF5 . 4 ) 


F5 . 3 ’ ) 


INITIALIZE THE BOUNDRY CONDITIONS 


VC1)=0. 

V C N ) =0 . 

N1=N-1 

DO 30 1=2, N1 
30 VCI)=100. 

CALL EXP FIN DIF FOR INFINITE HEAT TRANSFER COEFFICENT 


OQO OQO OQQ OOO OQO 


90 


CALL EFDIHC ( N , NTOT , TSI , V , T , KR , BETA ) 

STOP 

END 

SUBROUTINE EFDIHC 

SUBROUTINE EFDIHC ( N , NTOT , TSI , V , T , KR , BETA ) 

FOR INFINITE HEAT TRANSFER COEFFICENT 

IMPLICIT REAL*8 CA-H,0~Z) 

REAL*8 VTC100),VC100),M(100),P(100) 

REAL*8 KR 
N1=N-1 


BEGIN TIME STEP LOOP 
DO 20 L=1 , NTOT 

CALCULATE DEPENDENT VARIABLE 

DO 40 1=2, N1 

IM1=I-1 

IP1=I+1 

VT(I)=VCI)+TSI*( (1 . +BETA*VCI) )*( VCIP1 )+V(IMl )-2 .*V(I) ) ) 
DERSQR=DABS( V ( IP1 ) -V( IM1 ) ) 

IF ( DEP.SQR . LE . 0 . 0 ) GO TO 40 

VT ( I ) = VT ( I ) +TSI* ( ( BETA* ( DERSQR ) **2 . )/4. ) 

40 CONTINUE 

PRINT RESULTS 

WRITEC6 ,5) 

5 FORMATC/) 

WRITE ( 6 , 31 )L 

31 FORMATC IX, 'TIME STEP NUMBER=' ,13/) 

TIME=T*DFLOAT(L) 

WRITE (6, 32) TIME 

32 FORMAT (5X, 'ELAPSED TIME= ', FI 0 . 4 ,' SECONDS '/ ) 

DO 70 1=1, N 

70 V(I)=VT(I) 

WRITE (6,81)CV(I),I=1,N) 

81 FORMATC 1X,11(F9.5,2X)) 

20 CONTINUE 
RETURN 
END 
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SOURCE. IMPVAR 

WRITTEN BY R.F. HANDSCHUH 

***** PROGRAM #6 ***** 

THIS ROUTINE IS TO BE USED FOR COMPAP.ISION TO EXPONENTIAL 
FINITE DIFFERENCE ALGORITHM. THIS ROUTINE WILL USE THE 
IMPLICIT ROUTINE TO SOLVE FOR THE TEMPERATURE FIELD USING 
THE TRI-DIAGONAL MATRIX ALGORITHM. 


IMPLICIT REAL*8CA-H»0-Z) 
REAL*8 V ( 1 0 0 ) , KR 


INPUT PROGRAM DATA 
WRITE ( 6,15) 

15 FORMAT( IX, 'NUMBER OF NODES=N I3V) 

READ( 9, 10 )N 

10 FORMATCI3 ) 

WRITE (6,16) 

16 FORMAT C IX, 'TOTAL NUMBER OF TIME STEPS= NTOT 13') 

READ( 9,21 )NTOT 

21 FORMATCI3 ) 

WRITEC 6,29) 

29 FORMATC IX, 'INPUT TIME*THERMAL DIFFUSIVITY / LENGTH SQUARED 
READ( 9,25 )TSI 

25 FORMAT (F5. 3) 

WRITE(6 ,22) 

22 FORMATC IX, 'TOTAL TIME OF ONE TIME STEP= T F7 . 5 ' ) 

READC 6,23 )TS 

23 FORMATC F7 . 5 ) 

WRITE ( 6,19) 

19 FORMATC IX, 'INPUT REFERENCE THERMAL CONDUCTIVITY=F5 . 9 ' ) 
READC 9 , 17 )KR 

17 FORMAT (F5. 9) 

WRITEC 6,26) 

26 FORMATC IX, 'INPUT THERMAL CONDUCTIVITY SLOPE VALUE= F5.9’) 
READC 9,27 )BETA 

27 FORMATCF5 . 9 ) 


INITIALIZE BOUNDRY CONDITIONS 


VC1) = 0 . 

V(N)=0 . 

N1=N-1 

DO 30 1=2, N1 
30 VCI)=100. 

CALL IMPLICIT ROUTINE 


F5 . 3 ' ) 
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CALL IMPL ( N , NTOT , TSI , V , TS , KR , BETA ) 

STOP 

END 

SUBROUTINE IMPL 

SUBROUTINE IMPL ( N , NTOT , OMEGA , T , TS , KR , BETA ) 

FOR INFINITE HEAT TRANSFER COEFFICENT AND VARYING THERMAL CONDUCTIVITY 
IMPLICIT REAL*8 CA-H,0-Z) 

REAL*8 AC101),BC101),CC101),DC101),KAP(101),GAM(101),TC101) 

REAL*8 TOC 101) 

REAL*8 KR 

PRINT HEADING 

WRITEC 6,331) 

331 FORMATC1X, ****** SOURCE . IMPVAR *****•/) 

RHO= 1 . 0 
CP=1 . 0 
N1=N-1 

BEGIN TIME STEP LOOP 
DO 20 L=1 , NTOT 

CALCULATE THOMAS ALGORITHM VARIABLES AND THOSE THAT ARE A FUNCTION 
OF TEMPERATURE 

DO 21 1=1, N 
21 D(I)=T(I) 

DO 200 1=1, N 
200 TO(I)=T(I) 

DO 25 1=2, N1 
KAP(I)=1 . 0+BETA*TCI) 

25 GAM(I)=BETA*CTCI+1)-T(I-1 ))/(«. ) 

DO 30 1=2, N1 

A(I)=-OMEGA*(KAP(I)-GAM(I) ) 

B ( I ) = ( 1 . +2.*OMEGA*KAP(I) ) 

30 CCI)=-OMEGA*CKAPCI)+GAMCI) ) 

CALL TRI-DIAGONAL-MATRIX ALGORITHM 

CALL TRIDAGC 2 ,N1,A,B,C,D,T) 

0 CONTINUE 

PRINT THE RESULTS 

WRITEC 6, 5) 

5 FORMAT C / ) 


z'** 

L, 
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WRITE( 6 , 31 )L 

31 FORMAT ( IX, 'TIME STEP NUMBER 3 ' ,13/) 

TIME=TS*DFLOAT(L) 

WRITE( 6 , 32 )TIME 

32 F0RMAT(5X, 'ELAPSED TIME 3 FI 0 . A , ' SECONDS '/ ) 

WRITE ( 6,81)(T(I),I=1,N) 

81 FORMAT ( 1X,11(F9.5,2X)) 

20 CONTINUE 
RETURN 
END 

SUBROUTINE TRIDAG 

THIS ROUTINE IS FOR THE SOLUTION OF THE THOMAS ALGORITHM 

THIS ROUTINE WAS TAKEN FROM THE BOOK APPLIED NUMERICAL METHODS 
BY CARNAHAN, LUTHER, AND WILKES. 

SUBROUTINE TRIDAG (IF,L,A,B,C,D,V) 

IMPLICIT REAL*8 C A-H , O-Z ) 

REAL*8 A ( 1 0 1 ) , B ( 1 0 1 ) , C ( 1 0 1 ) , D ( 1 0 1 ) , V ( 1 0 1 ) , BETA (101), GAMMA (101) 

COMPUTE INTERMEDIATE ARRAYS BETA AND GAMMA 

BETA(IF)=B(IF) 

GAMMA ( IF )=D( IF) /BETA (IF) 

IFP1=IF+1 
DO 1 I=IFP1,L 

BETA(I)=B(I)-A(I)*C(I-1)/BETA(I-1) 

1 GAMMA ( I )=(D(I)-A( I )*GAMMA ( 1-1 ) )/BETA( I) 

COMPUTE FINAL SOLUTION VECTOR V 

V(L)=GAMMA(L) 

LAST=L-IF 
DO 2 K 3 1 , LAST 
I = L-K 

2 V(I)=GAMMA(I)-C(I)*V(I+1 )/BETA(I) 

RETURN 

END 
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SOURCE. COUE 

WRITTEN BY R.F. HANDSCHUH 
***** PROGRAM #7 ***** 


THIS PROGRAM IS TO BE USED TO DEMONSTRATE THE USE OF THE 
EXPONENTIAL FINITE DIFFERENCE METHOD ON THE DEVELOPING TEMPERATURE 
FIELD IN A LAMINAR COUETTE FLOW 


IMPLICIT REAL*8(A-H,0-Z) 
REAL*8 V(IOO) 

INPUT PROGRAM DATA 


WRITEC 6,15) 

15 FORMAT( IX, ’NUMBER OF NODES=N 13’ /) 

READ ( 9 , 1 0 )N 

10 FORMATCI3 ) 

WRITEC 6,12) 

12 FORMAT ( IX, 'NUMBER OF SUB INTERVALS= NS 13*) 

P.EADC 6,13 )NS 

13 FORMAT (13) 

WRITEC 6,16) 

16 FORMATC IX, ’TOTAL NUMBER OF POSITION STEPS= NTOT 13') 

READ! 9,21 )NTOT 

21 FORMATC 13 ) 

WRITEC 6, 24) 

24 FORMATC IX, 'INPUT (POSITION STEP*KIN. VISO/CLENGTH SQUARED) F5.3’) 
READC 9 , 25 )TSI 

25 FORMATCF5 . 3) 

WRITEC 6,22) 

22 FORMATC IX, 'TOTAL DISTANCE OF ONE STEP= T F5.3’) 

READC 6 , 23 ) T 

23 FORMATCF5 . 3) 

WRITEC 6,14) 

14 FORMATC IX, 'INPUT (MAX WIDTH)/(FLOW VEL) F5.4’) 

READ ( 6 , 1 7 )SP 

17 FORMAT (F5. 4) 

WRITEC 6,31) 

31 FORMATC IX, ’INPUT INTERVAL FOR PRINTING RESULTS= 13') 

READC 6 , 32 )IPR 

32 FORMATC 13) 


INITIALIZE THE BOUNDRY CONDITIONS 

VC 1 ) = 0 . 

DO 30 1=2, N 
V(I)=100 . 


30 
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CALL EXPONENTIAL FINITE DIFFERENCE FOR COUETTE FLOW 

CALL EFFL (N,NS, NTOT , TSI , V , T , SP , IPR ) 

STOP 

END 

SUBROUTINE EFFL 

SUBROUTINE EFFL ( N , NS , NTOT , TSI , V , T , SP , IPR ) 

FOR COUETTE LAMINAR FLOW 
IMPLICIT REAL*8(A-H,0-Z) 

REAL*8 VT(IOO) ,V(100) ,M( 100 ) ,P( 1 00 ) ,FL(20 ) 
TS=TSI/DFLOAT( NS+ 1 ) 

N1=N-1 


PRINT HEADING 
WRITE (6,92) 

92 FORMAT ( IX, ’SOLUTION FOR DEVELOPING TEMPERATURE FIELD IN’ 

* f LAMINAR COUETTE FLOW'//) 

DY=1 ./DFLOAT(N-l) 

CALCULATE PARAMETER THAT VARIES WITH POSITION IN THE FLOW 
DO 115 1=2, N 

115 FL ( I ) =SP/ ( DY*DFLOAT ( I- 1 ) ) 

NS1 =NS+1 

BEGIN TOTAL POSITION STEP LOOP 

DO 20 L=1 , NTOT 
IT=IT+1 

ZERO THE SUM OF DRIVE NUMBERS 

DO 15 1=1, N 
15 P(I)=0. 

SET TEMPORARY VALUES EQUAL TO THE LAST POSITION STEP VALUE 

DO 10 1=1, N 
10 VT( I ) =V( I ) 

SUB POSITION INTERVAL 

DO 30 K= 1 , NS1 

CALCULATE THE DRIVE NUMBERS 

DO 40 1=2, N1 
IM1=I-1 
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IP1=I+1 

4 0 M(I)=(2.*VTCI)-VTCIM1)-VTCIP1))/VTCI) 

CALCULATE TEMPORARY DEPENDENT VARIABLE ON SUB-INTERVAL 
DO 50 11=2, N1 

50 VT(I1 )=VT(I1 )*DEXPC -TS*MCI1 )*FL(I1 ) ) 

DO 60 1=2, N1 
60 P(I)=P(I)+M(I) 

30 CONTINUE 

CALCULATE THE COMPLETE STEP DEPENDENT VARIABLES 
DO 70 1=1, N 

V(I)=V(I) *DEXP( -TS#PCI)*FLCI)) 

70 CONTINUE 

IF( IT . LT . IPR )GO TO 20 
IT=0 

PRINT THE RESULTS 

WRITE( 6,5) 

5 FORMATC / ) 

WRITEC 6 , 31 )L 

31 FORMAT(lX, 'POSITON STEP NUMBERS, 13) 

DIST=T*DFLOAT ( L ) 

WRITEC 6 , 32 )DIST 

32 FORMATC5X, 'LOCATION DIST= ' , FI 0 . 4 , 'METERS ' ) 

WRITEC 6,82)(V(I),I=1,N) 

82 FORMAT C1X,11CF9.5,2X)) 

20 CONTINUE 
RETURN 
END 
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WRITTEN BY R. F. HANDSCHUH 
SOURCE. EX3D 

***** PROGRAM #8 ***** 


THIS ROUTINE IS FOR FINDING THE TEMPERATURE AT A GIVEN LOCATION 
AND TIME FOR A 3-DIMENSIONAL SOLID. 

THIS PROGRAM CALCULATES THE EXACT TEMP AS FOUND IN THE BOOK 
"TRANSPORT PHENOMENA" BY BIRD, STEWART, AND LIGHTFOOT . 


IMPLICIT P.EAL*8(A-H,0-Z) 

REAL*8 PI , PI2 , PI3 
PI=3 . 1415926 

INPUT PROGRAM DATA 

WRITE( 1,31) 

31 FORMAT ( IX, ’NUMBER OF NODES PER COORDINATE DIRECTION=I3 ' ) 
READ( 5,32 ) NODES 

32 FORMAT (13) 

WRITE( 1,1) 

1 FORMAT ( IX, 'INPUT INITIAL TEMPERATURE FOR THE SOLID F6.3') 
READC 5 , 2 )T0 

2 FORMAT(F6 . 3 ) 

WRITE (1,3) 

3 FORMAT( IX, ’SURFACE TEMPERATURE FOR TIME > 0 SECONDS F6.3') 
READ( 5 , 9 )T1 

9 FORMAT (F6. 3) 

WRITE ( 1,5) 

5 FORMAT( IX, 'INPUT THERMAL DIFFUSIVITY F5.3') 

READ (5,6) ALFA 

6 FORMAT (F5. 3) 

ASSUMING 3 DIM. CUBE 

DATA A, B, C/1. 0,1. 0,1.0/ 

WRITE( 1,9) 

9 FORMAT( IX, 'INPUT THE NUMBER OF TIMES THROUGH SUMMATION 13') 
READ( 5 , 1 1 )N1 

11 FORMAT(I3 ) 

WRITE ( 1,12) 

12 FORMAT( IX, 'INPUT TIME TO BE EVALUATED AT F5.3') 

READ( 5 , 1 9 )TIME 

19 FORMAT (F5. 3) 

WRITE ( 6,20)T0,T1 

20 FORMAT( IX, 'TEMP INITIAL= ' , El 5 . 8 , 2X , ' TEMP SURF= ' , E15 . 8/ ) 

WRITE (6,30) ALF A , TIME 

30 FORMAT(lX, 'THER DIFF= ’ , El 5 . 8 , 2X , ' TIME= ' , F7 . 5 ) 
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START SUMMATION ROUTINE 

PI3=PI**3. 

PI2=PI**2. 


CALCULATE THE TEMPERATURE ALONG THE DIAGONAL 

DEL=A/(DFL0AT(N0DES-1 ) ) 

DO 100 NODE=l, NODES 
X=DEL*DFLOAT(NODE~l) 

Y=X 

SUM=0 . 0 

DO 10 K=1,N1 

RP=DFLOAT(K-l) 

RP1 =RP+ . 5 

C0SP=DC05( RP1*PI*Z/C ) 

DO 10 N=1,N1 
RN=DFLOAT(N-l) 

RN1 =RN+ . 5 

COSN=DCOS(RNl*PI*Y/B) 

DO 10 M=1 , N1 
RM=DFLOAT(M-l) 

RM1 =RM+ . 5 

COSM=DCOS ( RM1*PI*X/A ) 

J=M+N+K-3 
Jl=J/2 
J2=J1*2 
VAL=-1 . 0 

IF( J.EQ. J2)VAL=1 . 0 

GAM= ( RM1**2 ) / ( A*A ) + ( RN1**2 )/ ( B*B ) + ( RP1**2 ) / ( C*C ) 
EXPLIM= ( -GAM*PI2* ALFA* TIME ) 

IF ( EXPLIM .LT.-100. )GO TO 999 
EP=DEXP( EXPLIM) 

GO TO 998 
999 CONTINUE 
EP= 0 . 0 

998 CONTINUE 

SUM=SUM+ C VAL/( RM1*RN1*RP1*PI3 ) ) *EP*COSM*COSN*COSP 
10 CONTINUE 

TNEW=T1+8.*(T0-T1)*SUM 

PRINT THE RESULTS 

WRITE (6,16)X,Y,Z 

16 FORMAT! IX, ' X= ’ , F5 . 3 , 2X , ’ Y= ' , F5 . 3 , 2X , ’ Z= ’ , F5 . 3 ) 

IFCTHEW.LT. 1E-10)TNEW=0 . 0 
WRITE (6,15 )TNEW 

15 FORMAT! IX, ’ TEMPER ATURE= ' , El 5 . 8 ) 

100 CONTINUE 
STOP 
END 
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SOURCE. EFD3D 

WRITTEN BY R.F. HANDSCHUH 
***** PROGRAM #9 ***** 


THIS PROGRAM IS FOR 3-DIMENSIONAL CARTESIAN COORDINATES 

UNSTEADY STATE HEAT TRANSFER IN A CUBE. THE METHOD OF SOLUTION IS THE 

EXPONENTIAL FINITE DIFFERENCE ALGORYTHM. THIS PARTICULAR PROGRAM IS 

FOR INFINITE HEAT TRANSFER COEFFICENT AT THE EXPOSED SURFACES 

AT X=Y=Z=1 . 0 FOR THE THREE SURFACES WHERE X,Y,Z EQUAL 0.0 

ARE TO BE CONSIDERED AS PERFECTLY INSULATED. 


IMPLICIT REAL*8 C A-H, O-Z) 

REAL*8 VC25,25,25) 

INPUT PROGRAM DATA 

WRITEC6 , 15) 

15 FORMATC IX, ’NUMBER OF NODES=N 13'/) 

READ ( 9 , 1 0 ) N 

10 FORMATCI3 ) 

WRITEC6, 12) 

12 FORMATC IX, 'NUMBER OF TIME SUB INTERVALS 3 NS 13') 

READ ( 9 , 1 3 ) NS 

13 FORMATCI3 ) 

WRITE ( 6 , 16) 

16 FORMATC IX, 'TOTAL NUMBER OF TIME STEPS 3 NTOT 13’) 

READ C 9 , 2 1 ) NTOT 

21 FORMATCI3 ) 

WRITEC 6,24) 

24 FORMATC IX, 'INPUT TIME*THERMAL DIFFUSIVITY / LENGTH SQUARED F5.3’) 
READC 9 , 25 )TSI 

25 FORMAT (F5. 3) 

WRITEC 6, 22) 

22 FORMATC IX, 'TOTAL TIME OF ONE TIME STEP 3 T F6 . 4 ' ) 

READC 9 , 23 )T 

23 FORMAT (F6. 4) 

WRITEC 6, 26) 

26 FORMATC IX, 'NUMBER OF TIME STEPS BEFORE PRINTING 3 13') 

READC 9 , 27 )IPR 

27 FORMATCI3 ) 

WRITEC 6 ,250 )N, NS, NTOT 

250 FORMAT C IX , ' # OF NODES 3 ', 13 , 2X, ’ f OF SUB-TIME-INT 3 ' ,13 , 2X, 

*'# OF TIME STEPS 3 ’,13) 

WRITEC6 , 251 )TSI,T 

251 FORMATC IX, ' (TIME*THER DIFF)/LENGTH SQUARED 3 ' F5 . 3 , 2X , 

* ’ TIME STEP LENGTH 3 ' ,F6. 4/) 

N1=N-1 


INITIALIZE BOUNDRY CONDITIONS 


no non non non non non non 
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DO 30 

1=1, N1 


DO 30 

J=1,N1 


DO 30 

K= 1 , N 1 

30 

V(I, J 
I = N 

,K)=1 .0 


DO 50 

J = 1,N 


DO 50 

K= 1 , N 

50 

V C I , J 
J=N 

, K ) = 0 . 0 


DO 51 

1=1 ,N 


DO 51 

K= 1 , N 

51 

V( I , J 
K = N 

, K ) = 0 . 0 


DO 52 

1 = 1, N 


DO 52 

J= 1 , N 

52 

V ( I , J 

, K ) = 0 . 0 


CALL EXP FIN DIF FOR INFINITE HEAT TRANSFER COEFFICENT 

CALL EFDIHC ( N , NS , NTOT , TSI , V , T , IPR ) 

STOP 

END 

SUBROUTINE EFDIHC 

SUBROUTINE EFDIHC ( N , NS , NTOT , TSI , V , T , IPR ) 

FOR INFINITE HEAT TRANSFER COEFFICENT 
IMPLICIT REAL*8 (A-H,0-Z) 

REAL*8 VT(25,25,25) ,V(25,25,25) , M( 25 , 25 , 25 ) , P ( 25 , 25 , 25 ) 
REAL*8 Ml, M2 
TS=TSI/DFLOAT ( NS+ 1 ) 

N1=N-1 
NS1 =NS+1 

PRINT HEADING 

WRITEC 6,200 ) 

200 FORMATC1X, *********** RESULTS FROM EFD3D **********•//) 
START TOTAL TIME STEP LOOP 
DO 20 L= 1 , NTOT 

ZERO THE SUM OF THE SUB-INTERVAL DRIVE NUMBERS 

DO 15 K=1 , N 
DO 15 J = 1 , N 
DO 15 1=1, N 
15 P ( I , J , K ) =0 . 

SET SUB-INTERVAL VALUES EQUAL TO THE LAST TIME STEP VALUES 


non non ooooo 


101 


DO 10 K=1,N 
DO 10 J=1,N 
DO 10 1=1, N 

10 VT(I,J,K)=V(I,J,K) 

SUB TIME INTERVAL 

CALCULATE THE DRIVE NUMBERS WHICH IS DEPENDENT ON LOCATION IN THE CUBE 

DO 30 KS= 1 , NS1 
DO 42 K=2,N1 
KM1 =K- 1 
KP 1 = K + 1 
DO 41 J=2 , N1 
JM1 = J- 1 
JP1=J+1 
DO 40 1=2, N1 
IM1 =1“ 1 
IP1=I+1 

Ml =VT ( IP1 , J , K ) + VT ( IM1 ,J , K ) +VT( I , JP1 , K ) + VT ( I , JM1 , K ) 

M2 =M 1 + VT C I , J , KP 1 ) + VT ( I , J , KM1 ) - 6 . *VT ( I , J , K ) 

IF(VT(I,J,K) . LE.0.0)M(I,J,K)=0. 

IF ( VT ( I , J , K ) . LE . 0 . 0 ) GO TO 40 
MCI, J,K)=M2/VT(I, J,K) 

40 CONTINUE 

41 CONTINUE 

42 CONTINUE 

INSULATED BOUNDRY ALONG X-AXIS 


J = 1 
K= 1 

KP1=K+1 
JP1=J+1 
DO 48 1=2, N1 
IP1=I+1 
IM1 =1- 1 

M1=VT(IP1, J,K)+VT(IM1, J,K)+2.*VT(I, JP1 , K ) +2 . *VT(I , J,KP1) 
IFCVTCI, J,K) .LE. 0 . 0 )M(I, J,K)=0 . 

IF( VT( I , J , K ) . LE . 0 . 0 )GO TO 48 
M(I,J,K)=(Ml-6 . *VT( I ,J,K))/VTCI,J,K) 

48 CONTINUE 

INSULATED BOUNDRY ALONG Y-AXIS 


1=1 

IP1=I+1 
K= 1 

KP1 =K+ 1 
DO 43 J=2,N1 
JP1=J+1 
JM1 = J- 1 

Ml =2 . *VT ( IP1 , J , K ) +VT ( I , JP1 , K ) +VT( I , JM1 , K ) +2 . *VT( I , J , KP1 ) 
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IFCVTCI, J,K) .LE. 0 . 0 )MCI, J,K)=Q . 

IFCVTCI, J,K).LE.O.O)GO TO 43 

MCI, J,K)=CMl-6 .*VTCI, J,K) )/VT(I, J,K) 

43 CONTINUE 

INSULATED BOUNDRY ALONG Z-AXIS 
J=1 

JP1=J+1 

1=1 

IP1=I+1 
DO 44 K=2 , N1 
KM1=K- 1 
KP1 =K+ 1 

M1=2.*VTCIP1 , J,K)+2.*VTCI, JP1 ,K)+VTCI, J,KP1)+VTCI, J,KM1) 

IFCVTCI, J,K) .LE. 0 . 0 )MCI, J,K)=0. 

IFCVTCI, J,K).LE.O.O)GO TO 44 

MCI, J,K)=(ni-6 .*VTCI, J,K) )/VTCI, J,K) 

4 CONTINUE 

INSULATED FACE AT Z=Q 

K= 1 

KP1=K+1 

DO 45 1=2, N1 

IP1=I+1 

IM1=I-1 

DO 45 J=2,N1 

JP1=J+1 

JM1 = J- 1 

M1=VTCIP1 ,J ,K) +VT C IM1 , J , K ) +VT C I, JP1 , K ) +VTC I , JM1 , K ) +2 . *VT C I , J , KP1 ) 
IFCVTCI, J,K) .LE. 0 . 0 )M(I, J,K)=0 . 

IFCVTCI, J,K).LE.O.O)GO TO 45 
MCI,J,K)=CMl-6. *VTCI, J ,K) )/VTCI, J,K) 

45 CONTINUE 

AT THE FACE WHERE Y=0 
J=1 

JP1=J+1 

DO 46 1=2, N1 

IP1=I+1 

IM1=I-1 

DO 46 K=2,N1 

KP1=K+1 

KM1=K-1 

M1=VTCIP1, J,K)+VTCIM1, J,K)+2.*VTCI, JP1,K)+VTCI, J,KP1 )+VTCI, J,KM1) 
IFCVTCI, J,K) .LE.O. Q)MC I, J,K)=0. 

IFCVTCI, J,K).LE.O.O)GO TO 46 

MCI, J,K)=CMl-6 .*VTCI, J,K))/VTCI, J,K) 

46 CONTINUE 

AT THE FACE WHERE X=0 
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1=1 

ipi=i+i 

DO 47 J=2,N1 

JP1=J+1 

JM1=J-1 

DO 47 K=2 , N1 

KP1=K+1 

Kni=K-l 

M1=2.*VT(IP1, J,K)+VT(I, JP1,K)+VT(IV JM1 , K ) +VTCI, J, KP1 ) +VTCI, J, KM1 ) 
IFCVTCI, J,K) .LE. 0 . 0)11(1, J,K)=0. 

IFCVTCI, J,K).LE. 0.0) GO TO 47 

MCI, J,K)=(Ml-6 .*VTCI, J,K) )/VT(I, J,K) 

47 CONTINUE 

CORNER AT ORIGIN 

Ml = 2 . *VT( 1 , 2 , 1 ) +2 . *VT( 2 , 1 , 1 ) +2 . *VT( 1 , 1 , 2 ) -6 . *VT( 1,1,1) 
M(1,1,1)=M1/VTC 1,1,1) 


CALCULATE THE SUB-INTERVAL DEPENDENT VARIABLES 

DO 50 11=1, N 
DO 50 Jl=l , N 
DO 50 K1=1,N 

IF(M(I1,J1,K1) .LT.-50. )VTCI1,J1,K1)=0.0 
. IFCMCIl , J1 ,K1 ) . LT . -50 )GO TO 50 

VTCI1. J1,K1)=VTCI1,J1,K1)*DEXP(TS*MCI1,J1,K1)) 

50 CONTINUE 

SUM THE DRIVE NUMBERS 

DO 60 1=1, N 
DO 60 J=1 , N 
DO 60 K= 1 , N 

60 PCI, J,K)=P(I,J,K)+MCI,J,K) 

30 CONTINUE 

CALCULATE THE NEXT COMPLETE TIME STEP DEPENDENT VARIABLES 

DO 70 K=1,N 
DO 70 J=1 , N 
DO 70 1=1, N 

IFCPCI,J»K) .LT.-50. )VCI,J,K)=0. 0 

IFCPCI, J,K) .LT.-50 . )GO TO 70 

VCI, J,K)=VCI, J,K)*DEXPCTS*PCI, J,K)) 

70 CONTINUE 

ITMAX=ITMAX+1 

IF C UMAX . LT . IPR ) GO TO 2 0 

WRITEC 6 , 5 ) 

5 FORMATC / ) 

WRITEC 6 , 31 )L 

31 FORMATC IX, 'TIME STEP NUMBER= ' ,13/ ) 


nnn 


I on 


TIt1E=T*DFL0AT C L ) 

WRITEC 6 » 32 )TIME 

32 FORMAT (5X, 'ELAPSED TIME= \ FI 0 . 4 SECONDS V ) 

PRINT OUT THE DIAGONAL RESULTS 

WRITEC 6 ,82)CV(I,I,I)>I=1>N) 

82 FORMATC 1 1 C 2X, F8 . 6 ) ) 

ITMAX= 0 
20 CONTINUE 
RETURN 
END 
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SOURCE. EXPL3D 

WRITTEN BY R.F. HANDSCHUH 

***** PROGRAM #10 ***** 


THIS PROGRAM IS FOR 3-DIMENSIONAL CARTESIAN COORDINATES 

UNSTEADY STATE HEAT TRANSFER IN A CUBE. THE METHOD OF SOLUTION IS THE 

PURE EXPLICIT METHOD. THIS PARTICULAR PROGRAM IS 

FOR INFINITE HEAT TRANSFER COEFFICENT AT THE EXPOSED SURFACES 

AT X=Y=Z=1 . 0 FOR THE THREE SURFACES WHERE X,Y,Z EQUAL 0.0 

IS TO BE CONSIDERED AS PERFECTLY INSULATED. 


IMPLICIT REAL*8(A-H,0-Z) 

REAL*8 V(25,25,25) 

INPUT PROGRAM DATA 

NUMBER OF NODES = N 

READ( 5 , 1 0 )N 
10 FORMAT (13) 

TOTAL NUMBER OF TIME STEPS = NTOT 

READ( 5 , 21 )NTOT 
21 FORMAT (14) 

TIME*THERMAL DIFFUSIVITY/LENGTH SQUARED = TSI 

READ( 5 , 25 )TSI 
25 FORMAT( F5 . 3 ) 

TOTAL TIME OF ONE TIME STEP = T 

READ( 5 , 23 )T 
23 FORMAT(F6 . A ) 

NUMBER OF STEPS BEFORE PRINTING THE RESULTS = IPR 

READ( 5 , 27 )IPR 
27 FORMAT( 13 ) 

WRITE (6,252) 

252 FORMAT( IX , ******************* PURE EXLICIT FINITE DIFFERENCE ’ 
* ’ MFTHOD ***************** • /// ) 

WRITE ( 6 , 250 )N, NS , NTOT 

250 FORMAT ( IX, '# OF NODES= ' , 13 , 2X , ' # OF SUB-TIME-INT= ' , 13 , 2X , 

*'# OF TIME STEPS=',I4) 

WRITE ( 6 , 25 1 )TSI , T 

251 FORMAT( IX, ' (TIME^THER DIFF)/LENGTH SQUARED='F5. 3,2X» 

* ' TIME STEP LENGTH= ' , F6 . A/ ) 

N1=N-1 
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INITIALIZE BOUNDRY CONDITIONS 

DO 30 1=1, N1 
DO 30 J=1,N1 
DO 30 K=1,N1 
30 V(I,J,K)=1.0 

I=N 

DO 50 J=1,N 
DO 50 K=1,N 

50 V(I,J,K)=0.0 
J = N 

DO 51 1=1, N 
DO 51 K=1,N 

51 VCI,J,K)=0.0 
K = N 

DO 52 1=1, N 
DO 52 J=1,N 

52 V(I,J,K)=0.0 

CALL PURE EXPLICIT FINITE DIFFERENCE FOR INFINITE HEAT TRANSFER COEFFICENT 

CALL PURE(N,NTOT,TSI,V,T,IPR) 

STOP 

END 

SUBROUTINE PURE 

SUBROUTINE PUREC N , NTOT , TSI , V , T , IPR ) 

PURE EXPLICIT FINITE DIFFERENCE METHOD 
FOR INFINITE HEAT TRANSFER COEFFICENT 

IMPLICIT REAL*8(A-H,0-Z) 

REAL*8 V(25,25,25) , VT(25,25,25) 

REAL*8 Ml, M2 
N1=N-1 

START TIME STEP LOOP 
DO 20 L= 1 , NTOT 

SAVE VALUES FROM THE LAST TIME STEP 

DO 39 1=1, N 
DO 39 J=1,N 
DO 39 K=1,N 

39 VT(I, J,K)=V(I, J,K) 

CALCULATE THE FIELD VARIABLE USING THE EXPLICIT FINITE DIFFERENCE METHOD 

DO 42 K=2,N1 
KM1=K-1 
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KP1=K+ 1 

DO 41 J=2,N1 

JM1 = J- 1 

JP1=J+1 

DO 40 1=2, HI 

IM1=I-1 

IP1=I+1 

M1=VT(IP1,J,K)+VT(IM1,J,K)+VT(I,JP1 ,K)+VT(I, JM1 ,K) 

M2=M1 + VT ( I , J , KP1 ) +VTC I , J , KM1 ) -6 . *VT ( I , J , K) 

VC I, J,K)=V(I, J,K)+TSI*M2 

40 CONTINUE 

41 CONTINUE 

42 CONTINUE 

INSULATED BOUNDRY ALONG X-AXIS 

J=1 

K=1 

KP1=K+1 

JP1=J+1 

DO 48 1=2, N1 

IP1=I+1 

IM1=I-1 

M1=VT(IP1, J,K)+VT(IM1, J,K)+2.*VT(I, JP1 ,K)+2 .*VT(I, J,KP1 ) 
M2=Ml-6 . *VT(I, J ,K) 

VCI, J,K)=V(I, J,K)+TSI*ri2 
48 CONTINUE 

INSULATED BOUNDRY ALONG Y-AXIS 
1=1 

IP1=I+1 

K=1 

KP1=K+1 
DO 43 J=2,N1 
JP1=J+1 
JM1= J- 1 

Ml =2 .*VT(IP1 , J,K)+VT(I, JP1,K)+VT(I, JM1 ,K)+2.*VT(I, J,KP1) 
M2=Ml-6 .*VT(I, J,K) 

V(I,J,K)=V(I,J,K) +TSI*M2 

43 CONTINUE 

INSULATED BOUNDRY ALONG Z-AXIS 


J=1 

JP1=J+1 

1=1 

IP1=I+1 
DO 44 K=2,N1 
KM1=K- 1 
KP1 =K+ 1 

Ml =2 .*VT(IP1, J,K)+2.*VT(I, JP1,K)+VT(I, J , KPl ) +VT( I , J,KM1) 
M2=Ml-6 .*VT(I, J,K) 

V(I,J,K)=V(I,J,K) +TSI*M2 
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4 4 CONTINUE 

INSULATED FACE AT Z=0 
K=1 

KP1=K+1 

DO 45 1=2, N1 

IP1=I+1 

IM1=I-1 

DO 45 J=2,N1 

JP1=J+1 

JM1= J- 1 

M1=VT(IP1 , J,K)+VT(IM1 , J,K)+VT(I, JP1 ,K)+VT(I, JM1 ,K)+2 .*VT(I, J,KP1 ) 
M2=t11-6 .*VT(I, J,K) 

V(I,J,K)=V(I,J,K) +TSI*M2 

45 CONTINUE 

AT THE FACE WHERE Y=0 
J=1 

JP1=J+1 

DO 46 1=2, N1 

IP1=I+1 

im=M 

DO 46 K=2 , N1 

KP1=K+1 

KM1 =K- 1 

M1=VT(IP1, J,K)+VT(It11,J,K)+2.*VT(I, JP1,K)+VT(I, J, KP1 ) +VTCI , J, KM1 ) 
M2=M1 -6 . *VT ( I , J , K ) 

V(I, J,K)=V(I, J,K)+TSI*M2 

46 CONTINUE 

AT THE FACE WHERE X=0 
1=1 

IP1=I+1 
DO 47 J=2,N1 
JP1=J+1 
JM1=J- 1 
DO 47 K=2,N1 
KP 1=K+ 1 
KM1=K-1 

M1=2.*VT(IP1, J,K)+VT(I, JP1,K)+VT(I,JM1,K)+VTCI, J,KP1)+VT(I, 
t12=Ml-6.*VT(I, J,K) 

VC I, J,K)=V(I, J,K)+TSI*M2 

47 CONTINUE 

CORNER AT ORIGIN 

M1=2.*VT(1,2,1)+2.*VT(2,1,1)+2.*VT(1,1,2)-6.*VT(1,1,1) 

V(1,1,1)=V(1,1,1)+M1*TSI 


ITMAX=ITMAX+ 1 
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IF(ITMAX.LT. IPR)G0 TO 20 
WRITEC 6,5) 

5 FORMATC/) 

WRITEC 6 , 31 )L 

31 f ORMATC IX, 'TIME STEP NUMBER 3 ' , 13/ ) 
TIME=T*DFLOAT ( L ) 

WRITEC 6 , 32 )TIME 

32 FORMATC 5X, ’ELAPSED TIME 3 ' , F 1 0 . 4 , ' SECONDS ' / ) 

PRINT OUT THE DIAGONAL RESULTS 

WRITEC 6, 82 )( VC I, I, I), 1=1, N) 

82 FORMATC 11C2X,F8.6) ) 

ITMAX=0 
20 CONTINUE 
RETURN 
END 
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SOURCE. DOUGLA 


WRITTEN BY R.F. HANDSCHUH 
***** PROGRAM #11 ***** 

THIS PROGRAM IS FOR 3-DIMENSIONAL CARTESIAN COORDINATES 

UNSTEADY STATE HEAT TRANSFER IN A CUBE. THE METHOD OF SOLUTION IS THE 

METHOD OF DOUGLAS. THIS PARTICULAR PROGRAM IS 

FOR INFINITE HEAT TRANSFER COEFFICENT AT THE EXPOSED SURFACES 

AT X=Y=Z=1 . 0 FOR THE THREE SURFACES WHERE X,Y,Z EQUAL 0.0 

IS TO BE CONSIDERED AS PERFECTLY INSULATED. 


IMPLICIT REAL*8( A-H.O-Z) 

P.EAL*8 T( 25 , 25 » 25 ) 

INF UT PROGRAM DATA 

WRITE ( 6 , 15 ) 

15 FORMAT C IX, ’NUMBER OF NODES=N 13’/) 

READ ( 5 > 1 0 )N 

10 FORMAT (13) 

WRITEC6, 16) 

16 . FORMAT ( IX, ’TOTAL NUMBER OF TIME STEPS= NTOT 13’) 

READ( 5, 21 )NTOT 

21 FORMAT (13) 

WRITE( 6,24) 

24 FORMAT( IX, ’INPUT TIME*THERMAL DIFFUSIVITY / LENGTH SQUARED F5 . 3 ' ) 
READ( 5,25 )TSI 

25 FORMAT ( F5 . 3 ) 

WRITE( 6,22) 

22 FORMAT( IX, 'TOTAL TIME OF ONE TIME STEP= T F6.4') 

READ (5,23 )DT 

23 FORMAT(F6 . 4 ) 

WRITE( 6,26) 

26 FORMAT ( IX, ’NUMBER OF TIME STEPS BEFORE PRINTING= 13’) 

READ( 5 , 27 )IPR 

27 FORMAT( 13 ) 

WRITE ( 6 , 2 5 0 ) N , NS , NTOT 

250 FORMATUX,'# OF NODES= ’ , 13 , 2X , ’ # OF SUB-TIME-INT= ' , 13 , 2X , 

*’# OF TIME STEPS= ’ ,13) 

WRITE(6 ,251 )TSI , DT 

251 FOPMAT( IX , ’ (TIME*THER DIFF)/LENGTH SQUARED 3 ’ F5 . 3 , 2X , 

* ’ TIME STEP LENGTH= ' , F6 . 4/ ) 

N1=N-1 


INITIALIZE BOUNDRY CONDITIONS 


DO 30 1=1, N1 
DO 30 J=1 , N1 
DO 30 K=1,N1 
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30 T(I,J,K)=1.0 

C 

I=M 

DO 50 J=1 , N 
DO 50 K = 1,N 

50 T(I,J,K)=0.0 
J=N 

DO 51 1=1, N 
DO 51 K=1,N 

51 TCI,J,K)=0.0 
K = N 

DO 52 1=1, N 
DO 52 J=l, N 

52 TCI , J , K ) =0 . 0 
C 

CALL DG ( N , NS , NTOT , TSI , T , DT , IPR ) 

stop 

END 

SUBROUTINE DG 

SUBROUTINE DG C N , NS , NTOT , TSI , T , DT , IPR ) 

FOR INFINITE HEAT TRANSFER COEFFICENT 
IMPLICIT REAL*8(A-H,0-Z) 

REAL-8 U(25,25,25) , VC 25, 25, 25) ,M(25,25,25) ,T(25,25,25) 

REAL*8 TEMP (25) ,A(25) ,BC25) ,C(25) ,D(25) 

REAL*8 Ml, M2 
N1=N-1 

PRINT OUT HEADING 
WRITE ( 6,200) 

200 FORMATC1X, *********** RESULTS FROM METHOD OF DOUGLAS *********** - 
8 //) 

CALCULATE COEFFICENTS FOR THOMAS ALGORITHM ( TRI-DIAGONAL MATRIX SOLVER) 

AC 1 ) = 0 . 0 
B( 1 ) =1 . +TSI 
C(1)=-TSI 
DO 60 1=2, N 
ACI)=-.5*TSI 
BCD = 1 . +TSI 
60 C ( I ) =A ( I ) 

BEGIN TIME STEP LOOP 

DO 20 L=1 , NTOT 

CALCULATE TEMPORARY VARIABLES "U" - X DIRECTION, f, "V" - Y DIRECTION, 

THEN ACTUAL FIELD VARIABLE "T" - Z DIRECTION. 
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DO 

120 

LOOP=l,3 

DO 

42 

K=2,N1 

DO 

41 

J=2,N1 

DO 

40 

1=2, N1 


IF( LOOP . EQ . 2 )G0 TO 140 
IF ( LOOP . EQ . 3 )GO TO 141 
DELX=T(I+1,J,K)+T(I-1,J,K)-2.*T(I,J,K) 

DELY=T (I,J+1,K)+T(I,J-1,K)-2.*T(I,J,K) 

DELZ=T(I, J,K+1)+T(I, J,K-1)-2.*T(I, J,K) 

M(I,J,K)=T(I,J,K) +TSI* ( . 5*DELX+DELY+DELZ ) 

GO TO 40 

140 CONTINUE 

DELU=U (I+1,J,K)+U(I-1,J,K)-2.*U(I,J,K) 
DELY=T(I,J+1,K)+T(I, J-1,K)-2.*T(I, J,K) 
M(I,J,K)=M(I,J,K)+. 5*TSI* ( DELU-DELY) 

GO TO 40 

141 CONTINUE 

DELV=V (I,J+1,K)+V(I,J-1,K)-2.*V(I,J,K) 

DELZ=T( I ,J,K+1)+T(I,J,K-1)-2.*T(I,J,K) 
M(I,J,K)=M(I,J,K)+. 5*TSI* C DELV-DELZ ) 

40 CONTINUE 

41 CONTINUE 

42 CONTINUE 

INSULATED BOUNDRY ALONG X-AXIS 

J=1 
K = 1 

KP1=K+1 
JP1=J+1 
DO 48 1=2, N1 
IP1=I+1 
IM1=I-1 

IF (LOOP . EQ . 2 )GO TO 148 
IF ( LOOP . EQ . 3 ) GO TO 248 

M1=.5*(T(I+1,J,K)+T(I-1, J,K)-2.*T(I, J,K) ) 

M2=2.*T(I, JP1,K)+2.*T(I, J , KP1 ) -4 . *T ( I , J,K) 
mi,J,K)=T(I,J,K)+TSI*(fU+M2) 

GO TO 48 
148 CONTINUE 

DELU = U (I+l,J,K)+U(I-l,vJ,K)-2.*U(I,J,K) 
DELY=2.*(T(I,J+1,K)-T(I, J,K) ) 

MCI, J,K)=tt(I, J,K) + . 5*TSI*(DELU-DELY) 

GO TO 48 
248 CONTINUE 

DELV=2 .*V(I,J+1,K)-2.*V(I,J,K) 

M(I, J,K)=M(I,J,K)+.5#TSI*(DELV-2.*(T(I,J,KP1)-T(I, J,K))) 
48 CONTINUE 

INSULATED BOUNDRY ALONG Y-AXIS 


1=1 

I.P1 = I+1 
K=1 


non ooo 


113 


KP1=K+1 
DO 4 3 J=2,N1 
JP1=J+1 
JM1=J-1 

IF ( LOOP . EQ . 2 )G0 TO 143 
IF (LOOP . EQ . 3 )G0 TO 243 

m=T(IPl, J,K)+T(I, JP1,K)+T(I, Jt11,K)+2.*T(I, J,KP1 )-5 . *T(I, J,K) 
MI, J,K)=T(I, J,K)+TSI*M1 
GO TO 43 

143 CONTINUE 

DELU=2.*(U(IP1, J,K)-U(I, J,K) ) 

DELY=T(I, JP1 ,K)+T(I,JM1,K)-2.*T(I,J,K) 

M( I , J , K ) =M( I , J , K ) + . 5*TSI* ( DELU-DELY ) 

GO TO 43 

243 CONTINUE 

DELV=V(I,JH,K)+V(I,J-1,K)-2.#V(I,J,K) 

MI, J,K)=MI, J,K)+. 5*TSI*(DELV~2 .*(T(I, J,KP1 )-T(I, J,K) ) ) 

43 CONTINUE 

INSULATED BOUNDRY ALONG Z-AXIS 
J = 1 

JP1=J+1 

1=1 

IP1=I+1 
DO 44 K=2 , N1 
KM1=K-1 
KP1=K+ 1 

IF( LOOP . EQ . 2 )GO TO 144 
IF ( LOOP . EQ . 3 ) GO TO 244 

M1=T(IP1, J,K)+2.*T(I, JP1,K)+T(I, J,KP1)+T(I, J,KM1 )-5 .*T(I, J,K) 
MI, J,K)=T(I, J,K)+TSI*M1 
GO TO 44 

144 CONTINUE 
M1=M(I, J,K) 

M2=M1+.5*TSI#(2.*(U(IP1, J,K)-U(I, J,K) )-2 .*(T(I,vJPl ,K)-T(I, J,K) )) 
MCI, J,K)=M2 
GO TO 44 

244 CONTINUE 

DELZ=T ( I , J , K+ 1 ) +T( I , J , K-l ) -2 . *T( I , J , K ) 

M( I , J , K ) =M( I , J , K) ♦ . 5*TSI* (2.*(V(I,JP1 ,K)-V(I,J,K) ) -DELZ ) 

44 CONTINUE 

INSULATED FACE AT Z=0 


K = 1 

KP1=K+ 1 
DO 45 1=2, N1 
IP1=I+1 

im=i-i 

DO 45 J=2,N1 

JP1=J+1 

JM1=J-1 

IF ( LOOP . EQ . 2 ) GO TO 145 
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IF ( LOOP . EQ . 3 ) GO TO 245 

M1=.5*CT(IP1, J,K)mini,J,K)-2.*T(I,J,K)) 

M2=T(I, JP1 ,K)+T(I, JM1 ,K)+2.*T(I, J , KP1 ) -4 . *T( I , J,K) 

H(I, J,K)=T(I, J,K)+TSI*(m+M2) 

GO TO 45 

145 CONTINUE 

DELU=U(I+1 ,J,K)+U(I-1,J,K)“2.*U(I,J,K) 

DELY=T(I, J+1,K)+T(I, J-l ,K)-2.*T(I, J,K) 
n(I,J,K)=n(I,J,K)+. 5*TSI* C DELU-DELY ) 

GO TO 45 

245 CONTINUE 

DELV=V(I,J+1,K)+V(I,J-1,K)-2.*V(I, J,K) 

M(I, J,K)=M(I, J,K)+.5*TSI*(DELV-2.*(TCI,J,KP1)-T(I, J,K) )) 

45 CONTINUE 

AT THE FACE WHERE Y=0 
J = 1 

JP1=J+1 

DO 46 1=2, N1 

IP1=I+1 

im=i-i 

DO 46 K=2 , N1 

KP1 =K+ 1 

KM1=K-1 

IF ( LOOP . EQ . 2 )GO TO 146 
IF ( LOOP . EQ . 3 ) GO TO 246 

ttl = .5*(T(IP1 , J,K)+T(It11, J,K)-2.*T(I, J,K) ) 
n2=2.*T(I>JPl,K)+T(I,J,KPl)+T(I,J,Kni)-4.*T(I,J,K) 

MCI, J,K)=T(I, J,K)+TSI*(M1+M2) 

GO TO 46 

146 CONTINUE 

DELU=U(I+1, J,K)+U(I-1, J,K)-2.*UCI,J,K) 

MI, J,K)=M(I, J,K)+.5*TSI*(DELU-2.*(T(I,JP1,K)-T(I, J,K)) ) 
GO TO 46 

246 CONTINUE 

DELZ=T(I,J,K+1)+T(I,J,K-1)-2.*T(I,J,K) 

MI, J,K)=MI,J,K)+.5*TSI*(2.*(V(I, JP1,K)-V(I, J,K) )-DELZ) 

46 CONTINUE 

AT THE FACE WHERE X=0 


1=1 

IP1=I+1 

DO 47 J=2,N1 

JP1=J+1 

JM1=J-1 

DO 47 K=2,N1 

KP1 =K+ 1 

KM1 =K- 1 

IF ( LOOP . EQ . 2 )GO TO 147 
IF C LOOP . EQ . 3 ) GO TO 247 

M1=T(IP1, J,K)+T(I, JP1 ,K)+TCI, J,KP1)+T(I, J,KM1)+T(I, JM1,K) 
MI, J,K)=T(I, J,K)+TSI*(M1-5.*T(I, J,K)) 
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GO TO 47 
147 CONTINUE 

DELY=T(I, J+1,K)+T(I,J-1,K)-2.*T(I, J,K) 

M(I, J,K)=M(I, J,K)+.5*TSI*(2.*(U(IP1, J,K)-U(I,J,K))-DELY) 

GO TO 47 
247 CONTINUE 

DELV=V(I, J+1,K)+V(I, J-1,K)-2.#V(I, J,K) 

DELZ=T(I, J,K+1)+T(I, J,K-1)-2.*T(I, J,K) 

M(I,J,K)=M(I,J,K)+. 5*TSI* ( DELV-DELZ ) 

47 CONTINUE 

COP.NER AT ORIGIN 

IF ( LOOP . EQ . 2 ) GO TO 151 
IF(L00P.EQ. 3)G0 TO 251 

M1=T(2,1,1)+2.*T(1,2,1)+2.*T(1,1,2)-5.*T( 1,1,1) 

M(1,1,1)=T (1,1,1 )+TSI*Ml 
GO TO 51 
151 CONTINUE 

mi=mc 1 , 1 , 1 ) 

M2=M1+.5*TSI*(2.*(U(2,1,1)-U(1,1,1))-2.*(T(1,2,1)-T(1,1,1))) 
M(1,1,1)=M2 
GO TO 51 
251 CONTINUE 

Ml=n( 1,1,1) 

M2=M1+.5*TSI*(2.*(V(1,2,1)-V(1,1,1))-2.*(T(1,1,2)-T(1,1,1))) 
M( 1 , 1 , 1 )=P12 
51 CONTINUE 


IF ( LOOP . EQ . 2 ) GO TO 130 
IF( LOOP . EQ . 3 )GO TO 230 

DO 70 K=1 , N1 
DO 70 J=1,N1 
DO 30 1=1, N1 
30 D(I)=M(I, J >K) 

CALL TRI DIAGONAL MATRIX ALGORITHM 

CALL TRIDAG (1,N1,A,B,C,D, TEMP ) 

C 

DO 430 1=1, N1 
430 U(I, J,K)=TEMP(I) 

C 

70 CONTINUE 
GO TO 330 
130 CONTINUE 

DO 71 K= 1 , N1 
DO 71 1=1, N1 
DO 72 J=1,N1 
72 D( J ) =M( I , J , K ) 

C 

CALL TRIDAG (1,N1,A,B,C,D, TEMP ) 
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DO 472 J=1,N1 
472 V(I, J,K)=TEMPC J) 

C 

71 CONTINUE 
GO TO 330 
230 CONTINUE 

DO 73 1=1, N1 
DO 73 J=1 , N1 
DO 74 K=1 , Ml 
74 D(K) 3 MCI,J,K) 

C 

CALL TRIDAG( 1,N1,A,B,C,D, TEMP ) 

C 

DO 474 K=1,N1 
474 TCI, J,K)=TEMPCK) 

C 

73 CONTINUE 
330 CONTINUE 
120 CONTINUE 

ITMAX=ITMAX+1 
IFCITMAX.LT. IPR)GO TO 20 
WRITEC 6,5) 

5 FORMATC/) 

WRITEC 6 , 91 )L 

91 FORMATC IX, 'TIME STEP NUMBER 3 ' ,13/) 

TIME=DT*DFLOATCL) 

WRITEC 6 , 32 )TIME 

32 FOP.MATC 5X, 'ELAPSED TIME 3 ’ , FI 0 . 4 ,' SECONDS '/ ) 

PRINT OUT THE DIAGONAL RESULTS 

WRITE C6,82)CT(I,I,I),I=1,N) 

82 FORMATC 11C2X,F8.6)) 

ITMAX 3 0 
20 CONTINUE 
RETURN 
END 

SUBROUTINE TRIDAG 

SUBROUTINE TP.IDAG CIF,L,A,B,C,D,V) 

IMPLICIT REAL*8CA-H,0-Z) 

REAL*8 AC 25) ,BC25),CC25) ,DC25) , VC 25 ) , BETAC 25 ) , GAMMAC 25 ) 

COMPUTE INTERMEDIATE ARRAYS BETTA AND GAMMA 

BETA(IF)=BCIF) 

GAMMA C IF )=DC IF) /BETAC IF) 

IFP1=IF+1 
DO 1 I=IFP1 , L 

BETACI) 3 BCI)-ACI)*CCI-1)/BETACI-1) 

1 GAMMACI)=CDCI)-ACI)*GAMMACI-1))/BETACI) 

C 


no 
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COMPUTE FINAL SOLUTION VECTOR V 

V( L) =GAMMA( L) 

LAST=L-IF 
DO 2 K=1 , LAST 
I=L-K 

2 VC I ) =GAMMA(I ) - C(I)*V(I+1 )/BETA(I) 

RETURN 
END 
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WRITTEN BY R.F. HANDSCHUH 

SOURCE. BURGER 

***** PROGRAM #12 ***** 

THIS PROGRAM IS FOR THE SOLUTION OF BURGER'S EQUATION 
BY THE EXPONENTIAL FINITE DIFFERENCE METHOD. 


IMPLICIT REAL*8(A-H,0-Z) 

REAL*8 V(100) 

READ IN DATA TO BE USED IN THE SOLUTION 


WRITE ( 6 , 15) 

15 FORMAT( IX, 'NUMBER OF NODES=N 13'/) 

READ( 9 , 1 0 )N 

10 FORMAT(I3) 

WRITEC6, 12) 

12 FORMATC IX, 'NUMBER OF TIME SUB INTERVALS= NS 13’) 
READ ( 9 , 1 3 )NS 

13 FORMAT (13) 

WRITE ( 6 , 16 ) 

16 FORMATC IX, 'TOTAL NUMBER OF TIME STEPS= NTOT 13') 
READ (9,21 )NTOT 

21 FORMATC 13 ) 

WP.ITEC 6,24) 

2A FOP.MAT( IX, 'INPUT TIME / LENGTH SQUARED F5.3') 

READ( 9 , 25 )TSI 

25 FORMAT (F5. 3) 

WRITE ( 6 , 26 ) 

26 FORMAT ( IX, 'INPUT KINEMATIC VISCOSITY= F5 . 3 ’ ) 

READ( 9,27 )RNU 

27 FORMAT(F5 . 3 ) 

WRITE (6, 22) 

22 FORMAT ( IX, 'TOTAL TIME OF ONE TIME STEP= T F5 . 3 ' ) 
P.EADC 9 , 23 )T 

23 FORMATC F5 . 3 ) 

WRITEC 6 , 1 A ) 

l<t FORMATC IX, 'INPUT NUMBUR OF STEPS BETWEEN PRINTS=I3’) 
READC 9,17 )IPR 

17 FORMATC 13) 

DATA FOR INITIAL AND BOUNDRY CONDITIONS 


30 


VC 1 ) = 0 . 

VC N ) = 1 . 
N1=N-1 

DO 30 1=2, N1 
V(I) = 1 . 


CALL EXPONETIAL FINITE DIFFERENCE FOR BURGER'S EQUATION 
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C 

CALL BURG ( N , NS , NTOT , TSI , V , T , RNU , IPR ) 

STOP 

END 

C 

C SUBROUTINE BURG 
C 

SUBROUTINE BURG (N, NS, NTOT, TSI, V,T, RNU, IPR) 

C 

C 

IMPLICIT REAL*8(A-H,0-Z) 

REAL*8 VT(100),V(100),M(100),P(100), THE (100) 

TS=TSI/DFLOAT(NS+l ) 

DX = 1 . O/DFLOAT ( N- 1 ) 

N 1 =N~ 1 
NS1=NS+1 
WRITE (6, 900) 

900 FORMAT (/****** SOLUTION FOR BURGER EQUATION *****'/) 

C 

C TOTAL TIME STEP LOOP 

C 

DO 20 L=1 , NTOT 
ITMAX=ITMAX+ 1 
C 

C ZERO THE SUM OF DRIVE NUMBERS 

C 

DO 15 1=1, N 
15 P(I)=0. 

C 

C SET THE TEMPOARY FIELD VARIABLE EQUAL TO THE LAST TIME STEP VALUE 

C 

DO 10 1=1, N 
10 VT( I ) =V ( I ) 

C 

C SUB TIME INTERVAL 

C 

DO 30 K=1 , NS1 
C 

C CALCULATE THE SUB-INTERVAL DRIVE NUMBERS 

C 

DO 40 1=2, N1 

IM1=I-1 

IP1=I+1 

IF(VTCI) .LE. 0 . 0)GO TO 40 

M( I ) =- 0 . 5*DX*(1 .-VT(I) )*(VT(IP1 )-VT(IMl) )/VT(I) 

M(I)=M(I)+RNU*( VTCIP1 )+VT(IMl )-2.*VT(I) )/VT(I) 

40 CONTINUE 

C 

C CALCULATE THE SUB-INTERVAL DEPENDENT VARIABLES 

C 

DO 50 11=2, N1 
CHECK=TS*M( II ) 

IF ( CHECK. LE. -50 . )VT(I1 )=0 . 0 
IF ( CHECK . LE . -5 0 . ) GO TO 50 
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VT ( II ) =VT( II ) *DEXP C TS*M( II ) ) 
50 CONTINUE 

SUM THE DRIVE NUMBERS 

DO 60 1=2, N1 
60 P( I ) =P( I ) +MC I ) 

30 CONTINUE 


CALCULATE THE DEPENDENT VARIABLE AT THE NEXT COMPLETE STEP 

DO 70 1=1, N 
CHECK=TS*P( I ) 

IF( CHECK. LE. -50. )V(I)=0.0 
IF( CHECK . LE . -50 . )GO TO 70 
V(I) =V(I)*DEXP(TS*P(I) ) 

70 CONTINUE 

OUTPUT THE RESULTS 

IF C ITMAX . LT . IPR ) GO TO 20 

ITMAX=0 

WRITE (6, 5) 

5 FORMAT ( / ) 

WRITE ( 6 , 31 )L 

31 FOP.MATdX, ’TIME STEP NUMBER=*,I3) 

TIME=T*DFLOAT ( L ) 

WRITE ( 6 , 32 )TIME 

32 FORMATC5X, 'ELAPSED TIME= ', FI 0 . 4 ,' SECONDS ' ) 

ISTEP=(N-1)/10 

DO 110 1=1, N 
110 THE(I)=1 . 0-V(I) 

DO 80 J=1 , ISTEP 
IS=( J-l )*11+1 
IFIN= J*1 0 + 1 

WRITEC6 ,81 ) (THE(I) ,I=IS,IFIN) 

81 FORMAT ( 1X,11(F8.6,2X)) 

80 CONTINUE 
20 CONTINUE 
RETURN 
END 
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WRITTEN BY R.F. HANDSCHUH 

SOURCE. EXBURG 

***** PROGRAM #13 ***** 


THIS PROGRAM IS FOR THE SOLUTION OF BURGER’S EQUATION USING AN EXPLICIT 
TECHNIQUE. THE RESULTS WILL BE USED TO COMPARE TO THE EXPONENTIAL 
FINITE DIFFERENCE TECHNIQUE. 


IMPLICIT REAL*8 C A-H, O-Z) 

REAL*8 V(IOO) 

INPUT PROGRAM DATA 

WRITE(6, 15) 

15 FORMATClX, ’NUMBER OF NODES=N 13'/) 

READ( 9 . 10 )N 

10 FORMAT( 13 ) 

WRITE ( 6,16) 

16 FORMATC IX, 'TOTAL NUMBER OF TIME STEPS= NTOT 13') 
READ (9,21) NTOT 

21 FORMATCI3) 

WRITE ( 6,24) 

24 FORMATC IX, 'INPUT TIME / LENGTH SQUARED F5 . 3 ' ) 

READC 9 , 25 )TSI 

25 FORMATCF5 . 3 ) 

WRITEC 6,26) 

26 FORMATC IX, 'INPUT KINEMATIC VISCOSITY= F5 . 3 ' ) 

READC 9 , 27 )RNU 

27 FORMATC F5 . 3 ) 

WRITEC 6,22) 

22 FORMATC IX, 'TOTAL TIME OF ONE TIME STEP= T F5.3') 

READ C 9 , 23 )T 

23 FORMAT (F5. 3) 

WRITEC 6 , 14) 

14 FORMATC IX, 'INPUT NUMBER OF STEPS BETWEEN PRINTS=I3') 
READC 9,17 )IPR 

17 FORMATC 13) 


30 


INITIALIZE THE BOUNDRY CONDITIONS 

VC 1 ) = 1 . 

V(N) = 0 . 

N1=N-1 

DO 30 1=2, N1 
VCI) = 0 . 


CALL EXPLICIT FINITE DIFFERENCE SOLUTION FOR BURGER'S EQUATION 

CALL BURG C N , NTOT , TSI , V , T, RNU , IPR ) 

STOP 
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END 

SUBROUTINE BURG 

SUBROUTINE BURG ( N , NTOT , TSI , V , T , RNU , IPR ) 


IMPLICIT REAL*8(A-H,0-Z) 

REAL*8 VTC100),VC100),THEC100) 

PRINT HEADING 

WRITE C 6,999 ) 

999 FORMATC IX, ***** EXPLICIT BURGER S EQT SOLUTION ****’/) 
DX=1.0/DFLOAT(N-1) 

N1=N-1 

TIME STEP LOOP 

DO 20 L=1 , NTOT 
ITMAX=ITMAXH 

SAVE LAST TIME STEP VALUES 

DO 10 1=1, N 
10 VT(I)=V(I) 

EVALUATE EXPLICIT FINITE DIFFERENCE EQUATION 

DO 40 1=2, N1 

IM1=I-1 

IP1=I+1 

VCI)=VTCI)-VTCI)*T*CVTCIP1)-VTCIM1) )/C2.*DX) 

V(I)=V(I) +RNU*T*( VTC IP1 ) +VTCIM1 ) - 2 . *VTCI) ) / C DX*DX ) 
wO CONTINUE 

WRITE OUT THE RESULTS 

IFdTMAX.LT. IPR) GO TO 20 
ITMAX= 0 
WRITEC 6,5) 

5 FORMATC/) 

WRITEC 6 , 31 )L 

31 FORMATC IX, ’TIME STEP NUMBER=’,I3) 

TIME=T*DFLOAT CL) 

WRITEC 6 , 32 )TIME 

32 FORMATC 5X, 'ELAPSED TIME= ' , F 1 0 . 4 , ' SECONDS ' ) 

ISTEP= ( N- 1 )/ 1 0 

DO 110 1=1, N 
110 THECI)=VCI) 

DO 80 J= 1 , ISTEP 
IS=C J-l )*11+1 
IFIN=U*10+1 

WRITEC 6 ,81) CTHECI) ,I=IS,IFIN) 
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81 FORMAT ( 1X,11(F8.6,2X)) 

80 CONTINUE 
20 CONTINUE 
RETURN 
END 
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SOURCE. NONBOU 

WRITTEN BY R.F. HANDSCHUH 

***** PROGRAM #14 ***** 


THIS PROGRAM IS USED FOR THE SOLUTION OF THE BOUNDRY LAYER FLOW 
OVER A FLAT PLATE. THE DIRECTION OF FLOW IS IN THE X-DIRECTION 
WHICH IS USED AS THE MARCHING DIRECTION FOR THE EXPONENTIAL FINITE 
DIFFERENCE ALGORITHM. THE THERMAL AND VELOCITY BOUNDRY LAYERS 
CAN BE EXTRACTED FROM THE TEMPERATURE AND VELOCITY FIELDS FOUND. 


IMPLICIT REAL*8(A-H,0-Z) 

REAL*8 U(101),T(101),V(101) 

INPUT PROGRAM DATA 

WRITEC6, 15) 

15 FORMAT( IX, 'NUMBER OF NODES IN Y DIRC=N 13'/) 

READC 9, 10 )N 

10 FORMATCI3) 

WRITEC6 , 12) 

12 FORMATC IX, 'NUMBER OF SUB INTERVALS 3 NS 13’) 

READC 9 , 1 3 )NS 

13 FORMAT (13) 

WRITE ( 6 , 16 ) 

16 FORMATC IX, 'TOTAL NUMBER OF X-DIR STEPS 3 NTOT 13') 
READC 9,21 )NTOT 

21 FORMAT (13) 

WRITEC 6,24) 

24 FORMATC IX, 'INPUT STEP LENGTH F5 . 3 ’ ) 

READC 9 , 25 )DX 

25 FORMATC F5. 3) 

WRITEC 6, 26) 

26 FORMATC IX, 'NUMBER OF STEPS BEFORE PRINTING 3 13') 
READC 9 , 27 )IPR 

27 FORMATCI3 ) 

WP.ITEC 6,110) 

110 FORMATC IX, 'INPUT KINEMATIC VISCOSITY 3 F6 . 4 ' ) 

READC 9,111) RNU 

111 FORMATCF6 . 4 ) 

WRITEC 6,101) 

101 FORMATC IX, 'INPUT THERMAL DIFFUSIVITY F6 . 4 ' ) 

READC 9,102 )RAL 

102 FORMATCF6 . 4 ) 

WRITEC 6,103) 

103 FORMATC IX, 'INPUT YMAX F5 . 1 ' ) 

READC 9, 104) YMAX 

104 FORMATCF5 . 1 ) 

WP.ITEC 6, 250 )N, NS, NTOT 

250 FORMATC IX, '# OF NODES 3 ', 13 , 2X # OF SUB-INT 3 ' , 13 , 2X, 
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*'# OF TIME STEPS=',I3) 

RN=DFL0AT(N-1 ) 

DY=YMAX/RN 
WRITE C 6 , 251 )DX,DY 

251 FORMATC IX, ’ DX= ' , F8 . 4 , 2X , ’ DY= ' , F8 . 4/ ) 

WRITE ( 6,252) RNU , RAL 

252 FORMATC IX, 'KINEMATIC VISCOSITY= ' , F6 . , * CM*CM/S ’ , 2X, 
# ' THERMAL DIFFUSIVITY= ' , F6 . 4 , ' CM*CM/S ' ) 

INITIALIZE BOUNDRY CONDITIONS 

DO 30 J=2 , N 
U( J) = l . 0 
V(«J) = 0.0 
30 T ( J ) = 1 . 0 

U(l)=Q.O 
T( 1 ) =0 . 0 


CALL NON 1 ( N , NS , NTOT , RNU , U , V , T , IPR , DX , DY , RAL 1 

STOP 

END 


SUBROUTINE NON1 C N , NS , NTOT , RNU , U, V,T,IPR,DX,DY,RAL) 


IMPLICIT REAL*8( A-H,0-Z) 

REAL*8 U ( 1 0 1 ) , MU (101),V(101),T(101), MT (101) 

REAL*8 PU ( 1 0 1 ) , PT ( 1 0 1 ) , UT ( 1 0 1 ) , TT ( 1 0 1 ) , VTC 1 0 1 ) 

REAL*8 THEC 101, 1000, 3), UT 1(101) 

PRINT HEADING 

WRITEC 6,5) 

WRITEC 6,222) 

222 FORMATC IX, ' #*##**##### SOURCE . NONBOU ************ ' // ) 

WRITEC 6, 223) 

223 FORMATC/, IX, ’SOLUTION FOR BOUNDRY LAYER FLOW PAST A FLAT PLATE’/) 
DY2=DY*DY 

TS=DX*RAL/ ( DY2*DFLOAT (NS+1 ) ) 

TSl=DX*RNU/(DY2*DFLOAT(NS+l ) ) 

DEL=DX/DFLOAT(NS+l) 

N1=N-1 

NS1 =NS+ 1 

NSTEP=(N-1)/10 

BEGIN TOTAL INTERVAL LOOP FOR L=1 TO NTOT STEPS 
DO 20 L= 1 » NTOT 

ZERO THE SUM OF DRIVE NUMBERS FOR THE NEXT SET OF SUB-POSITION INTERVALS 
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C 

DO 15 1=1, N 
PUCI)=0.0 
15 PTC I ) = 0 . 0 

C 

C SAVE LAST POSITION STEP VALUES FOR TEMPORARY VARIABLE CALCULATIONS 
C ON SUB-POSITION INTERVAL 

C 

DO 10 1=1, N 

VTC I ) =V ( I ) 

UT ( I ) =U ( I ) 

10 TT(I)=T(I) 

C 

C SUB - POSITION INTERVAL 

C 

DO 30 K=1,NS1 
C 

C CALCULATE TEMPERATURE FIELD DRIVE NUMBER 

C 

DO 41 J=2,N1 

JM1=J-1 

JP1=J+1 

MTCJ)=-VTCJ)*CTTCJP1)-TTCJM1) )*DY/C 2 . *RAL*UTC J)*TTC J) ) 

MTC J)=MT(J)+CTTCJP1)+TT(JM1)-2.*TT( J) )/CUTC J)*TTC J) ) 

41 CONTINUE 

C 

C CALCULATE X - DIRECTION VELOCITY DRIVE NUMBER 
C 

DO 141 J=2 , N1 
JM1 = J- 1 
JP1=J+1 

MUCJ)=-.5*VTCJ)*DY*CUTC JP1 )-UTC JM1 ) )/CRNU*UTC J)*UTC J) ) 

MUC J)=MUC J)+CUTCJP1)+UTC JM1)-2.*UTC J) )/(UT( J)*UTCJ) ) 

141 CONTINUE 

C 

C CALCULATE TEMPERATURE, X-DIRECTION VELOCITY, AND Y-DIRECTION VELOCITY 
C ON THE SUB-POSITION INTERVAL 

C 

DO 50 11=2, N1 

50 TT(I1)=TTCI1) *DEXP ( TS*MT C II ) ) 

C 

DO 51 1=2, N1 
UT1(I)=UT(I) 

51 UT(I)=UTCI)*DEXP(TS1*MU(I) ) 

C 

DO 65 J=2,N1 
JM1 = J-l 

VTC J ) =VTC JM1 ) - . 5*CDY/DEL)*CUTC J)-UT1 (J)+UTCJM1)-UT1(JM1)) 

65 CONTINUE 

C 

C SUM THE DRIVE NUMBERS 

C 

DO 60 J=2,N1 
PUC J ) =PU C J ) +MU C J) 
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60 PT ( J ) =PT C J ) +MT C J ) 

30 CONTINUE 

CALCULATE THE NEXT TOTAL POSITION STEP VALUES OF VELOCITIES AND TEMPERATURE 

DO 70 J=1,N 
UT(J)=U( J) 

UC J)=U( J)*DEXPCTS1*PUC J) )' 

T( J ) =T ( J ) *DEXPC TS*PTC J ) ) 

70 CONTINUE 

DO 75 J=2,N1 

75 VC J)=V( JM1)-.5*(DY/DX)*CUC J)-UTCJ)+UC JM1)-UTCJM1) ) 

ITMAX=ITMAX+1 

SAVE THE VALUES FOUND IN 3-DIMENSIONAL ARRAY "THE" 

DO 76 J=1,N 
THECJ,L,1)=UCJ) 

THE(J,L,2)=V(J) 

76 THE(J,L, 3)=T( J) 

IFCITMAX.LT. IPR)GO TO 20 

WRITE OUT THE RESULTS AT THE REQUESTED INTERVAL OF POSITION 

WRITE ( 6 , 5 ) 

5 FORMATC/) 

WRITEC 6 , 31 )L 

31 FORMAT (5X, 'POSITION STEP NUMBERS, 13) 

TSTEP=DX*DFLOATCL) 

WRITEC 6 , 32 )TSTEP 

32 FORMATC 5X , ’ X-POSITION= ' , FI 0 . 4/ ) 

WRITEC 6,101) 

101 FORMATC IX, 'THE U VELOCITY COMPONENT') 

DO 300 KK=1 , NSTEP 

IS=(KK-1 )*11+1 
IFIN=KK*1 0+KK 

WRITEC 6,82) CTHE(I,L, 1 ) ,I=IS,IFIN) 

300 CONTINUE 

82 FORMATC 11C2X,F8 .5) ) 

WRITEC 6,102) 

102 FORMATC IX, 'THE V VELOCITY COMPONENT') 

DO 301 KK=1, NSTEP 

IS=(KK-1 )*11+1 
IFIN=KK*1 0 +KK 

WRITE ( 6 , 82 ) C THE C I , L , 2 ) , I=IS , IFIN ) 

301 CONTINUE 
WRITEC 6,103) 

103 FORMATC IX, 'THE T FIELD VARIABLE ') 

DO 302 KK=1, NSTEP 

IS=(KK-1 )*11 + 1 
IFIN=KK*1 0 +KK 
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WRITE ( 6 , 8 2 ) ( THE C I , L , 3 ) , I=IS , IFIN ) 
302 CONTINUE 
ITMAX= 0 
20 CONTINUE 
RETURN 
END 
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